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ABSTRACT 

We  report  the  detection  of  far-IR  CO  rotational  emission  from  the  prototypical  Seyfert  2  galaxy  NGC  1068.  Using 
Herschel-PA.CS ,  we  have  detected  11  transitions  in  the  /uppcr  =  14-30  kp  =  580-2565  K)  range,  all 

of  which  are  consistent  with  arising  from  within  the  central  10"  (700  pc).  The  detected  transitions  are  modeled 
as  arising  from  two  different  components:  a  moderate-excitation  (ME)  component  close  to  the  galaxy  systemic 
velocity  and  a  high-excitation  (HE)  component  that  is  blueshifted  by  ~80  km  s_1.  We  employ  a  large  velocity 
gradient  model  and  derive  «h2  ~  1056  cm  ~3,  7km  ~  170  K,  and  Mh2  ~  106,7  M0  for  the  ME  component  and 
«H2  ~  106,4  cm-3,  7km  ~  570  K,  and  Mm  ~  105  6  M0  for  the  HE  component,  although  for  both  components 
the  uncertainties  in  the  density  and  mass  are  ±(0.6-0. 9)  dex.  Both  components  arise  from  denser  and  possibly 
warmer  gas  than  traced  by  low-/  CO  transitions,  and  the  ME  component  likely  makes  a  significant  contribution 
to  the  mass  budget  in  the  nuclear  region.  We  compare  the  CO  line  profiles  with  those  of  other  molecular  tracers 
observed  at  higher  spatial  and  spectral  resolution  and  find  that  the  ME  transitions  are  consistent  with  these  lines 
arising  in  the  ~200  pc  diameter  ring  of  material  traced  by  H2  1-0  S(l)  observations.  The  blueshift  of  the  HE  lines 
may  also  be  consistent  with  the  bluest  regions  of  this  H2  ring,  but  a  better  kinematic  match  is  found  with  a  clump 
of  infalling  gas  ~40  pc  north  of  the  active  galactic  nucleus  (AGN).  We  consider  potential  heating  mechanisms  and 
conclude  that  X-ray-  or  shock  heating  of  both  components  is  viable,  while  far-UV  heating  is  unlikely.  We  discuss 
the  prospects  of  placing  the  HE  component  near  the  AGN  and  conclude  that  while  the  moderate  thermal  pressure 
precludes  an  association  with  the  ~  1  pc  radius  H2O  maser  disk,  the  HE  component  could  potentially  be  located 
only  a  few  parsecs  more  distant  from  the  AGN  and  might  then  provide  the  Ay  ~  1025  cirT  2  column  obscuring 
the  nuclear  hard  X-rays.  Finally,  we  also  report  sensitive  upper  limits  extending  up  to  /upper  =  50,  which  place 
constraints  on  a  previous  model  prediction  for  the  CO  emission  from  the  X-ray  obscuring  torus. 
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1.  INTRODUCTION 

The  excited  molecular  gas  in  the  centers  of  Seyfert  galaxies 
offers  a  sensitive  probe  of  the  nature  of  active  galactic  nu¬ 
cleus  (AGN)  feedback  on  the  surrounding  interstellar  medium 
(ISM).  Observational  studies  of  the  most  highly  excited  ma¬ 
terial  in  Seyfert  nuclei  have  typically  used  the  H2  rotational 
(Lutz  et  al.  2000;  Rigopoulou  et  al.  2002;  Roussel  et  al.  2007) 
and  the  well-studied  rovibrational  (e.g.,  Thompson  et  al.  1978; 
Mouri  1994;  Maloney  1997;  Davies  et  al.  2005;  Rodrfguez- 
Ardila  et  al.  2005)  transitions.  The  pure  H2  rotational  lines 
( ±upper /^b  >  500  K)  are  easily  thermalized  at  moderate 
(n H2  >  103  cm~3)  densities,  while  the  rovibrational  lines 
(Fupper/^B  >  7000  K)  may  be  excited  through  collisions  in 
dense  («h2  A)  105  cm"  3)  gas  or  through  UV  fluorescence 
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(Sternberg  &  Dalgarno  1989).  Observations  of  these  tracers  in 
Seyferts  have  identified  a  number  of  potentially  important  exci¬ 
tation  mechanisms,  including  X-rays  from  the  AGN  (Maloney 
1997;  Rodriguez- Ardila  et  al.  2005);  shocks  associated  with 
supernova  (SN)  remnants,  radio  jets,  and  gravitational  instabil¬ 
ities  (Roussel  et  al.  2007);  and  stellar  far-UV  (FUV)  radiation 
(Davies  et  al.  2005),  with  no  clear  consensus  on  a  single  dom¬ 
inant  excitation  source.  The  far-IR  (FIR)  CO  rotational  tran¬ 
sitions  (CO [ /Upper  ->  /upper  ~  1],  with  /upper  «  13-50)  arise 
from  states  500-7000  K  above  ground,  have  critical  densities  of 
~106-108  cm  3,  and  complement  the  H2  transitions  for  stud¬ 
ies  of  warm  and  dense  material.  Compared  with  H2,  the  FIR 
CO  lines  trace  similar  energy  levels  but  have  higher  critical 
densities  and  are  less  sensitive  to  extinction.  Additionally,  the 
smaller  energy  gaps  between  levels  lead  to  a  finer  sampling  of 
density-temperature  phase  space.  These  lines  have  been  pro¬ 
posed  as  potential  tools  for  studying  the  obscuring  medium  of 
type  2  AGNs  (Krolik  &  Lepp  1989),  determining  the  energy 
budgets  of  composite  starburst/AGN  systems  (Meijerink  et  al. 
2007),  and  identifying  accreting  black  holes  in  the  early  uni¬ 
verse  (Spaans  &  Meijerink  2008;  Schleicher  et  al.  2010),  but 
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previous  facilities  were  unable  to  detect  this  line  emission  from 
extragalactic  sources.  Here,  we  take  advantage  of  the  superb 
sensitivity  of  Herschel-PACS  to  conduct  the  first  extragalac¬ 
tic  study  of  FIR  CO  emission,  from  the  prototypical  Seyfert  2 
galaxy  NGC  1068. 

NGC  1068  is  one  of  the  brightest  and  best-studied  Seyfert 
2  galaxies.  The  paradigm  of  an  optically  and  geometrically 
thick  molecular  torus  accounting  for  the  Seyfert  type  1  and 
2  dichotomy  followed  the  detection  of  scattered  broad-line 
emission  from  this  source  (Miller  &  Antonucci  1983),  and 
NGC  1068  has  been  at  the  center  of  subsequent  studies  of  the 
ISM  in  Seyfert  nuclei.  The  molecular  gas  in  the  central  V  of 
NGC  1068  has  been  well  studied,  and  here  we  review  some 
of  the  key  results.  Interferometric  observations  of  CO(l-O) 
have  identified  a  pair  of  ««15"  1 . 1  kpc)  radius  spiral  arms 

(Planesas  et  al.  1991;  Heifer  &  Blitz  1995;  Schinnerer  et  al. 
2000),  which  may  be  modeled  as  forming  in  response  to  a 
^17  kpc  bar  (Schinnerer  et  al.  2000).  These  arms  are  bright 
in  Br y  (Davies  et  al.  1998),  polycyclic  aromatic  hydrocarbon 
emission  (Le  Floc’h  et  al.  2001),  and  submillimeter  continuum 
(Papadopoulos  &  Seaquist  1999a)  and  contain  most  of  the  star 
formation  in  the  central  region.  Centered  on  the  AGN  is  the  ~5" 
(~350  pc)  circumnuclear  disk  (CND),  which  is  visible  in  CO 
and  H2  1-0  5(1)  but  becomes  particularly  prominent  in  images 
of  HCN  (Tacconi  et  al.  1994)  and  other  high-density  tracers 
(Garcfa-Burillo  et  al.  2010).  Strong  emission  in  CO(4-3)  and 
HCN(l-O)  indicates  that  the  gas  in  the  CND  is  both  warm  and 
dense  (Tacconi  et  al.  1994;  Sternberg  et  al.  1994;  Israel  2009; 
although  see  Krips  et  al.  201 1  for  a  lower  density  model).  The 
high  abundances  of  HCN,  CN,  H30+,  and  other  molecules  in 
the  CND  suggest  an  X-ray-driven  chemistry  (Usero  et  al.  2004; 
Garcia-Burillo  et  al.  2010;  Aalto  et  al.  2011),  and  X-ray  heating 
has  also  been  invoked  to  explain  the  strong  H2  1-0  5(1)  and 
[Fe  n]  emission  (Rotaciuc  et  al.  1991 ;  Maloney  1997;  Galliano  & 
Alloin  2002).  At  ~0’.'3  resolution,  the  H2  1-0  5(1)  observations 
resolve  the  CND  into  a  ~1"  ring-like  structure  (Galliano  & 
Alloin  2002),  while  the  line  spectral  profiles  show  evidence  for 
rotation,  expansion,  and  more  complex  kinematics  (Galliano  & 
Alloin  2002;  Galliano  et  al.  2003;  Davies  et  al.  2008).  Shocks 
following  this  non-circular  motion,  and  possibly  associated  with 
jet-ISM  interactions,  may  also  be  important  in  heating  the 
molecular  CND  (Krips  et  al.  2011).  At  ~0'.'l  resolution  the 
H2  1-0  5(1)  images  reveal  two  clumps  of  infalling  molecular 
material  at  ~0"  1  -iVA  scales  that  likely  play  an  important  role 
in  both  fueling  and  obscuring  the  AGN,  with  an  estimated  infall 
rate  of  ~15  M0  yr_1  to  within  a  few  parsecs  of  the  nucleus 
(Muller  Sanchez  et  al.  2009).  Finally,  milliarcsecond  resolution 
radio  observations  identify  a  series  of  H20  maser  spots  that 
trace  out  the  inner  surface  of  a  ~0.65  pc  radius  molecular  disk, 
centered  on  the  AGN  (Gallimore  et  al.  2004,  and  references 
therein). 

Hard  X-ray  observations  of  NGC  1068  indicate  large 
obscuration  to  the  nucleus  (Iwasawa  et  al.  1997;  Matt  et  al. 
1997;  Colbert  et  al.  2002)  by  an  intervening  medium  with  a  col¬ 
umn  density  possibly  exceeding  Ah  >  1025  cm-2  (Matt  et  al. 
1997).  Interferometric  mid-IR  observations  of  NGC  1068  iden¬ 
tify  a  parsec-scale  structure  of  hot  dust  that,  along  with  the  H20 
maser  disk,  may  represent  the  dusty  molecular  torus  responsible 
for  the  X-ray  obscuration  (Jaffe  et  al.  2004;  Raban  et  al.  2009). 
However,  other  investigators  have  found  evidence  that  at  least 
some  of  the  nuclear  obscuration  occurs  on  few  to  ten  parsec 
scales  from  the  AGN  (Cameron  et  al.  1993;  Honig  et  al.  2006; 
Muller  Sanchez  et  al.  2009). 


The  organization  of  this  paper  is  as  follows.  In  Section  2 
we  describe  the  Herschel-PACS  observations  of  the  FIR  CO 
lines  in  NGC  1068.  In  Section  3  we  analyze  the  gas  excitation 
and  estimate  physical  parameters,  and  in  Section  4  we  compare 
the  physical  parameters  and  line  profiles  with  those  of  other 
molecular  gas  tracers.  In  Sections  5  and  6  we  discuss  potential 
heating  mechanisms.  In  Section  7  we  discuss  our  detections  and 
upper  limits  in  the  context  of  the  molecular  ISM  within  a  few 
parsecs  of  the  AGN,  and  in  Section  8  we  summarize  our  findings. 
Throughout  this  paper,  we  adopt  a  distance  to  NGC  1068  of 
14.4  Mpc  (Bland-Hawthorn  et  al.  1997)  and  a  systemic  velocity 
Vlsr  =  1125  km  s'1. 

2.  OBSERVATIONS 
2.1.  Data  Acquisition  and  Reduction 

The  observations  were  made  with  the  Photodetector  Array 
Camera  and  Spectrometer  (PACS;  Poglitsch  et  al.  2010)  on 
board  the  Herschel  Space  Observatory  (Pilbratt  et  al.  2010). 
Most  of  the  data  presented  here  were  obtained  as  part  of  the 
SHINING  guaranteed  time  key  program.  The  SHINING  ob¬ 
servations  consisted  of  10  high-resolution  range  scans  concate¬ 
nated  to  cover  the  52-98  /xm  and  104-196  /xm  ranges,  as  well  as 
deeper  integrations  of  CO(17-16),  CO(24-23),  and  CO(40-39). 
These  latter  observations  targeted  CO  transitions  falling  in  rel¬ 
atively  clean  spectral  regions  and  were  conducted  to  provide 
a  coarse  but  sensitive  sampling  of  the  CO  spectral  energy  dis¬ 
tribution  (SED)  over  the  full  FIR  range.  The  SHINING  data 
yielded  detections  of  most  transitions  at  Juvvcr  V  24,  and  we 
obtained  follow-up  observations  of  CO(28-27)  and  CO(30-29) 
in  an  open  time  project  to  extend  our  CO  SED  measurements 
to  higher  J.  These  observations  amounted  to  a  total  of  13.7  hr 
of  integration  time.  The  data  reduction  was  done  using  the  stan¬ 
dard  PACS  reduction  and  calibration  pipeline  (ipipe)  included 
in  HIPE  5.0  975.  However,  for  the  final  calibration  we  normal¬ 
ized  the  spectra  to  the  telescope  flux  and  recalibrated  it  with  a 
reference  telescope  spectrum  obtained  from  dedicated  Neptune 
continuum  observations.  With  this  approach,  we  estimate  an 
absolute  flux  calibration  accuracy  of  30%. 

2.2.  Line  Flux  Estimation 

The  PACS  spectrometer  performs  integral  field  spectroscopy 
over  a  47"  x  47"  field  of  view,  resolved  into  a  5  x  5  array  of 
9'.' 4  spatial  pixels  (spaxels).  The  spectrometer  resolving  power 
varies  from  R  =  1000  to  3000  for  the  first-  and  second- 
order  observations  utilized  here.  In  Figure  1,  we  show  the 
spectra  from  the  central  spaxel  centered  on  11  of  the  12  CO 
transitions  falling  in  the  104-196  /xm  range.  The  CO(25-24) 
line  at  krest  =  104.44 /xm  lies  in  a  noisy  region  at  the  edge  of 
this  range  and  is  not  included.  Most  of  these  lines  are  strong  in 
the  central  spaxel,  but  little  flux  is  detected  outside  of  the  central 
spaxel,  as  expected  for  an  unresolved  source  {0SOUTCe  <  9'.'4).  All 
fluxes  and  upper  limits  presented  here  were  therefore  extracted 
from  the  central  spaxel  and  referenced  to  a  point  source  by 
dividing  by  the  recommended  point-source  correction  factors9 
(Poglitsch  et  al.  2010). 

The  CO  line  fluxes  were  measured  by  fitting  the  spectra  with  a 
Gaussian  profile  plus  a  baseline.  In  most  cases  a  linear  baseline 
was  adopted  and  the  three  parameters  defining  the  Gaussian 
were  allowed  to  vary  freely,  but  some  lines  required  a  modified 


9  See  also  http://herschel.esac.esa.int/twiki/pub/Public/PacsCalibrationWeb/ 
PacsSpectroscopyPerformanceAndCalibration_v2_4.pdf. 
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Figure  1.  Continuum-subtracted  spectra  of  the  7Upper  =  14—24  and  7Upper  =  30 
CO  lines.  The  blue  curves  show  the  line  fits,  and  features  other  than  CO  are 
labeled  in  red.  For  CO(16-15)  and  CO(17-16),  the  green  and  red  curves  show 
the  decomposition  of  the  fit  into  CO  and  other  lines,  respectively.  All  lines  are 
detected  with  the  exception  of  CO(23-22),  which  is  blended  with  a  strong  H2O 
line  209  km  s-1  to  the  red.  Here,  the  overplotted  green  curve  is  an  average  of 
the  CO(22-21)  and  CO(24— 23)  profiles  as  a  reference. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 
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Figure  2.  Top:  FIR  CO  line  fluxes  and  upper  limits  measured  here,  along 
with  lower-7  lines  from  the  literature.  The  brightest  set  of  7Upper  =  1—4  line 
fluxes  is  measured  in  11"-21"  beams  that  contain  a  mixture  of  the  CND  and 
the  more  extended  emission  (Israel  2009).  The  second  set  of  /upper  =  1-3 
points  is  interferometric  measurements  integrated  over  the  central  4"  from 
Krips  et  al.  (2011),  and  the  fainter  CO(l-O)  flux  is  the  same  from  Schinnerer 
et  al.  (2000).  The  shaded  region  shows  the  range  of  good-fitting  LVG  models 
(/r2d  =  X2/dof  1.1;  Section  3.2.3),  and  the  solid  black  curve  is  the 
single  best-fitting  model,  with  the  red  and  blue  curves  showing  the  individual 
contributions  from  the  ME  and  HE  components.  Middle:  central  velocities  of 
the  FIR  CO  lines,  with  average  values  of  the  ME  and  HE  line  centers  indicated 
with  horizontal  red  and  blue  lines.  Bottom:  measured  FWHM  of  the  FIR  CO 
lines,  with  tracks  indicating  the  expected  line  widths  (estimated  by  adding  the 
intrinsic  line  widths  and  the  spectrometer  resolution  in  quadrature )  for  sources 
with  intrinsic  widths  of  0,  100,  200,  and  400  km  s  1 . 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


approach.  A  broad  feature  underneath  the  CO(15-14)  line  is 
present  in  the  raw  data,  likely  due  to  an  imperfect  subtraction  of 
the  telescope  background,  and  we  remove  this  feature  using  a 
higher  order  baseline  fit.  The  integrated  CO(  15-14)  flux  and  the 
residual  line  profile  shown  in  Figure  1  are  consistent  with  those 
of  adjacent  transitions,  and  we  estimate  that  the  flux  uncertainty 
introduced  by  this  baseline  feature  is  less  than  the  assumed 
30%  absolute  flux  calibration  error.  The  CO(  16-15)  line  is 
blended  with  the  163  /xm  OH  doublet,  and  CO(17-16)  with 
a  pair  of  flanking  OH+  lines.  In  both  cases  we  estimate  the  CO 
flux  by  simultaneously  fitting  all  features.  An  unconstrained 
Gaussian  fit  to  the  relatively  low  signal-to-noise  ratio  (S/N) 
CO(22-21)  line  yields  a  much  broader  profile  than  for  any 
other  transition,  so  here  we  fix  the  width  of  the  CO  profile 
to  the  typical  value  derived  from  other  line  fits  (corresponding 
to  an  intrinsic  FWHM  of  250  km  s-1;  see  Figures  2  and  7  and 
the  discussion  in  Section  4.2.1).  CO(23-22)  is  blended  with  a 
strong  H2O  line  with  a  rest  wavelength  209  km  s_1  to  the  red. 
If  the  combined  feature  were  attributed  solely  to  H2O,  it  would 
be  both  broader  and  more  blueshifted  than  any  of  the  other 
six  H2O  lines  detected  in  the  PACS  scans,  and  we  interpret 
this  as  evidence  for  significant  contamination  by  CO(23-22).  A 
comparison  with  the  average  of  the  CO(22-21)  and  CO(24-23) 
profiles  suggests  that  the  CO(23-22)  line  is  weaker  and/or 
more  redshifted  than  these  lines  (Figure  1).  However,  due  to 
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Table  1 

PACS  CO  Line  Observations 


Line 

^•rest 

Om) 

Fluxa 

(10^17  Wm“2) 

Vob 

(km  s  J) 

FWHM 

(km  s-1) 

CO(14-13) 

186.00 

7.2  ±  2.3 

19  ±23 

305  ±  32 

CO(  15-14) 

173.63 

6.4  ±  2.2 

2  ±25 

313  ±37 

CO(16-15) 

162.81 

8.1  ±2.5 

49  ±25 

369  ±  28 

CO(17-16) 

153.27 

5.8  ±  1.8 

0  ±  26 

359  ±  27 

CO(18-17) 

144.78 

5.1  ±  1.6 

-15  ±28 

379  ±31 

CO(19-18) 

137.20 

2.6  ±0.9 

-31  ±35 

324  ±  55 

CO(20-19) 

130.37 

2.5  ±  0.9 

-62  ±  36 

297  ±  55 

CO(21-20) 

124.19 

2.4  ±0.9 

— 17  ±  46 

407  ±  87 

CO(22-21) 

118.58 

4.0  ±  1.4 

-44  ±  43 

387c 

CO(23-22) 

113.46 

Blended 

CO(24-23) 

108.76 

2.6  ±  1.0 

-94  ±  49 

385  ±  95 

CO(25-24) 

104.44 

<11.2 

CO(26-25) 

100.46 

n/a 

CO(27-26) 

96.77 

<5.8 

CO(28-27) 

93.35 

<4.6 

CO(29~28) 

90.16 

<9.3 

CO(30-29) 

87.19 

4.2  ±  1.9 

-89  ±50 

341  ±  117 

CO(31-30) 

84.41 

Blended 

CO(32-31) 

81.81 

<7.4 

CO(33~32) 

79.36 

<9.5 

CO(34-33) 

77.06 

<5.8 

CO(35-34) 

74.89 

<6.2 

CO(36-35) 

72.84 

<8.3 

CO(37-36) 

70.91 

<7.8 

CO(38~37) 

69.07 

<10.1 

CO(39-38) 

67.34 

<19.1 

CO(40-39) 

65.69 

<22.3 

CO(41^tO) 

64.12 

<13.2 

CO(42~41) 

62.62 

<21.6 

CO(43-42) 

61.20 

<16.8 

CO(44^13) 

59.84 

<10.8 

CO(45-44) 

58.55 

<14.8 

CO(46^15) 

57.31 

Blended 

CO(47^t6) 

56.12 

<14.7 

CO(48^17) 

54.99 

<28.3 

CO(49^t8) 

53.90 

<31.7 

CO(5CM9) 

52.85 

<45.2 

Notes. 

a  Total  uncertainties  combine  a  30%  calibration  error  with  statistical  errors  in 
line  fits.  Upper  limits  are  3 a  and  refer  to  the  flux  density  integrated  over  a 
600  km  s-1  bin.  Some  lines  are  blended  with  a  strong  feature,  and  it  is  not 
possible  to  obtain  a  flux  measurement  or  useful  upper  limit.  CO(26-25)  was  not 
covered  in  the  PACS  scans. 

b  Relative  to  Vlsr  =  1125  km  s-1.  Total  uncertainties  combine  a  spectral 
calibration  accuracy  of  10%  of  the  spectral  resolution  (20-30  km  s-1)  with 
statistical  errors  in  line  fits. 

c  Fixed  to  an  intrinsic  FWHM  of  250  km  s_1  (see  the  text). 


the  uncertainties  involved  in  deconvolving  the  CO  and  hCO 
lines,  we  simply  exclude  CO(23-22)  from  our  analysis.  All 
other  transitions  from  CO(14-13)  through  CO(24-23)  are  well 
detected  (Table  1). 

The  52-98  /tm  range  includes  the  CO(27-26)  through 
CO(50-49)  transitions.  None  of  these  lines  are  detected  in  the 
full  range  scans  or  the  targeted  CO(40-39)  observation  obtained 
with  SHINING,  but  our  follow-up  open  time  program  yielded  a 
detection  of  CO(30-29)  (Figure  1 ).  The  flux  for  this  line  was  es¬ 
timated  using  the  same  fitting  procedure  as  described  above  for 
the  ./upper  <  24  lines  (Table  1).  We  estimate  upper  limits  for  the 
non-detected  transitions  by  first  binning  the  data  to  600  km  s”1 
bins  and  then  estimating  the  3a  noise  levels  (Table  1). 


We  detect  no  emission  from  13CO.  The  13CO(  14-13)  transi¬ 
tion  at  A.rest  =  194.55  \i m  lies  at  the  noisy  edge  of  a  scan,  while 
the  /upper  >  2 1  transitions  at  A.  ^  1 04  p.  m  are  blended  with  1 2CO 
lines.  For  13CO(15-14)  through  13CO(20-19),  we  estimate  3a 
upper  limits  of  (2-4)  x  10”17  W  m  2,  for  a  600  km  s-1  bin 
size.  Our  most  stringent  lower  limit  on  the  12CO/13CO  flux  ra¬ 
tio  comes  from  the  /upper  =16  transition,  for  which  we  estimate 
I2CO(16— 15)/I3CO(  16— 15)  >  2.6.  Assuming  that  the  12COand 
13CO  SEDs  are  similar,  this  suggests  that  we  can  exclude  sig¬ 
nificant  contamination  of  the  detected  I2CO  lines  at  /upper  >21 
from  13CO. 

3.  EXCITATION  ANALYSIS 
3.1.  Evidence  for  Two  Components 

In  the  top  panel  of  Figure  2,  we  show  the  line  fluxes  and 
upper  limits  measured  here,  along  with  lower-/  measurements 
obtained  from  the  literature.  The  middle  and  bottom  panels 
show  the  central  velocities  and  widths  obtained  from  the 
Gaussian  fitting.  The  inflection  point  seen  in  the  FIR  CO 
line  SED  at  /upper  ~  19  suggests  the  presence  of  multiple 
components,  as  does  the  shift  in  central  velocities  between  the 
lowest-  and  highest-/  transitions.  For  simplicity,  we  assume  that 
the  FIR  CO  lines  are  produced  by  two  discrete  components: 
a  moderate-excitation  (ME)  component  near  the  systematic 
velocity  and  a  blueshifted  high-excitation  (HE)  component.  Our 
excitation  analysis  described  below  indicates  that  the  /upper  ^ 
17  and  /upper  Y  20  transitions  are  dominated  by  the  ME  and 
HE  components,  respectively.  Separately  averaging  the  central 
velocities  of  these  two  sets  of  lines  gives  Vme  =  17  ±  12  km  s'1 
and  Vhe  =  —  59  ±  20  km  s’  1  (Figure  2),  with  a  difference  of 
Vme  —  Ehe  =  76  ±  23  km  s_1 . 

3.2.  LVG  Modeling 

To  quantitatively  analyze  the  FIR  CO  line  SED,  we  employ 
a  large  velocity  gradient  (LVG)  model.  We  use  the  LVG 
calculation  described  in  Hailey-Dunsheath  et  al.  (2008),  with 
updated  CO-H2  collisional  coefficients  from  Yang  et  al.  (2010), 
and  a  thermalized  Hi  ortho/para  ratio.  In  this  model  the  shape  of 
the  CO  line  SED  is  determined  by  the  gas  density  («h2)>  kinetic 
temperature  (7kin),  and  CO  abundance  per  velocity  gradient 
([CO/Ht \/{dv/dr)).  The  source  is  assumed  to  consist  of  a 
number  of  unresolved  clouds,  and  the  absolute  line  luminosities 
scale  with  the  total  CO  mass  (Mco).  In  the  following  analysis, 
we  use  a  CO  abundance  of  [CO/H2]  =  I  O  ’  4  to  reparameterize 
[CO/H2]/(</u//r)  as  dv/dr  and  Mco  as  the  H2  mass  (MH 1). 
This  results  in  an  eight-parameter  model,  with  four  parameters 
for  each  of  the  ME  and  HE  components.  In  Section  3.2.5,  we 
discuss  the  effects  of  varying  the  CO  abundance. 

3.2.1.  Background  Radiation 

The  CO  excitation  may  be  affected  by  the  background 
radiation,  and  we  must  therefore  estimate  the  local  radiation 
density.  At  millimeter  wavelengths  the  background  is  dominated 
by  the  cosmic  microwave  background  (CMB;  Kamenetzky  et  al. 
2011),  but  the  FIR  background  arises  from  within  the  galaxy. 
The  PACS  integral  field  spectra  provide  a  continuum  map  of 
the  central  47"  x  47"  with  9'.'4  pixels  and  demonstrate  that  the 
continuum  emission  in  the  central  pixel  is  dominated  by  sources 
within  the  central  %  1 0" .  The  flux  density  in  the  central  pixel  may 
be  modeled  as  an  optically  thin  modified  blackbody  with  f  — 
1.5,  7jlisi  =  48  K,  and  normalized  to  FV(X  =100  p  m )  =  49  Jy, 
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and  we  adopt  this  as  an  estimate  of  the  continuum  brightness 
of  the  ~5"  CND.  If  the  flux  in  the  central  spaxel  is  indeed 
due  solely  to  the  CND,  then  this  is  a  moderate  (by  a  factor 
of  <2)  underestimate,  while  emission  from  within  the  central 
s»10"  but  outside  of  the  CND  may  also  be  contributing.  The 
measured  continuum  is  in  reasonable  agreement  with  previous 
calculations  and  observations.  For  comparison,  this  measured 
Fv  is  a  factor  of  I  -4  times  larger  in  the  X  =  50-200  / /  m  range 
than  obtained  from  the  radiative  transfer  modeling  in  Spinoglio 
et  al.  (2005).  Additionally,  extrapolating  our  modeled  SED  to 
X  =  450  /xm  yields  FV(X  =  450 /rm)  =1.2  Jy,  comparable  to 
the  peak  value  of  ~  1.5  Jy  beam-1  measured  in  a  ~9"  beam  by 
Papadopoulos  &  Seaquist  (1999a).  The  60  jrm/100  // m  ratio  in 
this  model  is  1.27,  and  the  FIR  flux10  is  2.7  x  10~12  W  m-2, 
giving  iFIR  =  4jrd2FFlR  =  1.7  x  101()  LQ. 

We  include  the  effects  of  background  radiation  on  the  equa¬ 
tions  of  statistical  equilibrium  following  Poelman  &  Spaans 
(2005).  The  important  parameter  in  this  approach  is  the  mean 
specific  intensity  of  the  external  radiation  field  at  the  cloud  sur¬ 
face,  which  we  define  as  Jv  ex t.  We  estimate  /y  ex t  using  a  simple 
geometrical  model  in  which  the  gas  clouds  are  uniformly  dis¬ 
tributed  in  a  sphere  with  an  observed  angular  size  Q.  and  are 
evenly  mixed  with  the  FIR-emitting  dust  grains.  For  optically 
thin  continuum,  the  mean  value  of  7„,ext  is  then  related  to  the 
observed  continuum  flux  density  Fv  o bS  as 


'v,ext 


=  iv. 


CB 


9  /yobs 

16  n  ’ 


(i) 


where  /„,cb  is  the  sum  of  the  CMB  and  cosmic  IR  background. 
We  take  £2  to  correspond  to  a  circular  diameter  of  4",  approx¬ 
imately  matched  to  the  size  of  the  CND  (see  Section  4  and 
Figure  8).  We  have  run  calculations  both  including  and  ignor¬ 
ing  the  local  contributions  to  the  background  and  see  negligible 
difference  in  the  results.  In  part  this  is  due  to  the  fact  that  the 
background  radiation  temperatures  are  only  7) ab  =  1 3-25  K  at 
the  wavelengths  of  the  detected  FIR  lines,  while  the  typical  exci¬ 
tation  temperatures  for  the  best-fitting  models  are  7’cx  ~  1 00  K 
and  Tex  ~  500  K  for  the  ME  and  HE  transitions,  respectively. 
In  addition,  most  of  the  lines  are  optically  thick  for  the  best¬ 
fitting  models,  and  hence  the  CO  is  insulated  from  the  external 
radiation  field. 


3.2.2.  Parameter  Limits 

We  explore  two-component  fits  to  the  FIR  CO  emission  over 
a  large  volume  of  eight-dimensional  parameter  space,  applying 
physical  limits  to  the  model  parameters.  The  most  important 
prior  restrictions  are  placed  on  the  velocity  gradient.  For  self- 
gravitating  clouds  in  virial  equilibrium,  we  can  approximate 
(dv/dr\u  ~  10  km  s-1  pc-1  («H2/105  cm-3)1/2  (Goldsmith 
2001).  The  actual  velocity  gradient  may  be  larger  due  to  addi¬ 
tional  sources  of  gravitational  potential,  a  high-pressure  inter¬ 
cloud  medium,  or  non-virialized  motion  (Bryant  &  Scoville 
1996),  but  smaller  values  are  unlikely.  Defining  Kv b  as  the  ratio 
between  dv/dr  and  ( dv/dr\u  (Papadopoulos  et  al.  2007) 

dv/dr  /  nm  \~1/2 
™  10  km s-1  pc-1  V 105  cm-3/  ’  U 

we  restrict  parameter  space  to  Kvlr  f  1 .  The  largest  measured 
velocity  gradient  in  NGC  1068  is  in  the  HjO  maser  disk 


10  Ffir  =  1.26  x  10-14  [2.58  /60/Jy  +  /100/Jy]  W  m-2. 


Table  2 

LVG  Model  Restrictions 


(1)  tfvir  ^  1 

(2)  dv/dr  1000  km  s-1  pc-1 

(3)  1.36  x  [MH2(ME)  +  MH2(HE)]  <  9  x  10s  Mq 

(4)  H2  rotational  lines  not  overproduced 


associated  with  the  AGN.  The  line-of-sight  velocities  of  the 
maser  spots  shift  by  ~600  km  s_1  over  a  ~2  pc  linear 
range  (Gallimore  et  al.  2001),  corresponding  to  an  effective 
dv/dr  ~  300  km  s_1  pc-1.  To  accommodate  the  maser  disk 
and  other  high-dispersion  regions  in  our  models,  we  extend  our 
calculations  up  to  dv/dr  =  1000  km  s_1  pc-1.  We  note  that 
restricting  Kv b  A  1  and  dv/dr  f  1000  km  s  1  pc-1  combine 
to  limit  the  density  to  hH2  4  109  cm-3,  but  as  we  discuss  below, 
such  high  densities  are  ruled  out  by  other  considerations.  We 
calculate  the  total  gas  mass  in  our  models  as  Mgas  =  1 .36  x  MH2, 
including  the  contribution  from  helium.  Schinnerer  et  al.  (2000) 
estimate  a  dynamical  mass  of  Mdyn  =  9x  108  Me  for  the  CND, 
and  we  discard  any  of  our  models  in  which  the  total  Mgas  of  the 
two  components  exceeds  Mdyn. 

The  range  of  parameter  space  allowed  by  the  CO  data  can 
be  further  reduced  by  considering  the  H2  pure  rotational  lines, 
which  arise  from  states  with  similar  upper  energy  levels  as  the 
FIR  CO  transitions.  We  calculate  the  H2  rotational  spectrum 
for  each  model,  under  the  simplifying  assumption  that  the  lines 
are  optically  thin  and  thermalized,  and  the  H2  ortho/para  ratio 
is  thermalized.  We  then  rule  out  any  model  that  overpredicts 
the  flux  in  any  of  the  lines  measured  in  the  large  (14"-27") 
apertures  of  the  Short  Wavelength  Spectrometer  (SWS)  onboard 
the  Infrared  Space  Observatory  (ISO)  (Lutz  et  al.  2000).  These 
prior  constraints  on  the  LVG  model  parameters  are  summarized 
in  Table  2. 

3.2.3.  General  Features  of  Good-fitting  Models 

We  proceed  by  generating  model  SEDs  over  a  regular  eight¬ 
dimensional  grid  and  for  each  model  calculating  /2  in  the 
normal  manner.  With  11  data  points  and  8  free  parameters, 
our  modeling  has  3  degrees  of  freedom  (dof).  Here,  we  dis¬ 
cuss  the  general  properties  of  the  set  of  solutions  for  which 
X2  —  Xmin  ^  1’  corresponding  to  x2/dof  4  1.1.  In  Figure  2,  we 
show  the  range  of  SEDs  covered  by  this  set  of  good  solutions, 
and  for  the  single  best-fit  model  we  show  the  decomposition  into 
the  ME  and  HE  components.  The  CO(18-17)  and  CO(19-18) 
lines  typically  receive  comparable  contributions  from  the  ME 
and  HE  components,  while  the  lower-  and  higher-7  transitions 
are  dominated  by  the  ME  and  HE  components,  respectively.  The 
shape  of  the  ME  SED  is  relatively  well  constrained  and  peaks 
in  the  7upPer  =  13-16  range.  The  single-dish  measurements  of 
CO(l-O),  CO(2-l),  CO(3-2),  and  CO(4-3)  we  show  in  Figure  2 
were  obtained  with  11"-21"  beams,  in  some  cases  comparable 
to  the  Herschel-PA.CS  resolution  (Israel  2009).  However,  these 
low-7  lines  receive  strong  contributions  from  the  lower  excita¬ 
tion  gas  in  the  »H5"  radius  starburst  ring  and  do  not  constrain 
our  models.  The  interferometric  measurements  of  the  CO(l-O) 
flux  in  the  CND  range  from  20  to  120  Jy  km  s_1  (Schinnerer 
et  al.  2000;  Krips  et  al.  2011),  larger  than  the  median  ME  model 
flux  of  7  Jy  km  s_1,  suggesting  that  the  ME  component  con¬ 
tributes  no  more  than  a  minor  fraction  of  the  observed  low-7 
emission.  The  HE  component  is  more  poorly  constrained  at  the 
high-7  end  and  is  well  fit  by  SEDs  peaking  from  7upper  =  25 
up  to  7upper  =  35.  For  the  latter  models,  our  upper  limits  to 
CO(28-27)  and  CO(34-33)  become  useful  constraints. 
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The  FIR  CO  emission  is  an  important  coolant  of  the  nuclear 
molecular  ISM  but  does  not  dominate.  The  total  luminosity 
emitted  in  the  11  transitions  detected  here  is  Leo  fir  =  3.3  x 
106lg,  and  summing  the  modeled  ME  and  HE  emission  over 
all  transitions  yields  Leo, me  +  Lco.he  =  (5.7-10.2)  x  106  L0. 
For  the  highest  excitation  models  some  additional  cooling  may 
arise  from  the  /upper  >  40  transitions  not  included  in  our  LVG 
calculation,  and  significant  emission  is  also  expected  in  the 
•Aipper  ^13  submillimeter  transitions.  The  total  emission  in  the 
H2  0-0  5(1),  5(3),  5(4),  5(5),  and  5(7)  rotational  lines  detected 
by  ISO-SWS  is  LUl  =  1 .7  x  107  L0  (Lutz  et  al.  2000).  Treating 
the  upper  limits  to  5(0),  5(2),  and  5(9)  as  detections  increases 
this  by  a  factor  of  1.5,  while  at  the  same  time  some  fraction 
of  the  lowest-/  emission  measured  with  the  largest  apertures 
may  arise  from  the  starburst  ring.  Our  PACS  scans  have  also 
detected  a  number  of  OH  and  H2O  transitions.  The  bulk  of  the 
emission  in  these  molecules  detected  in  the  central  spaxel  likely 
arises  from  the  unresolved  CND,  and  with  this  assumption  we 
estimate  nuclear  luminosities  of  Loh.fir  ^  1.5  x  107  L0  and 
Lh2o,fir  ^  3.0  x  106  L0.  The  FIR  range  includes  the  strongest 
OH  lines  at  79  /xm,  119/xm,  and  163  /xm,  while  for  H2O  (as 
with  CO)  the  longer  wavelength  emission  should  be  strong.  In 
the  FIR  the  CO  and  H2O  cooling  is  therefore  comparable,  while 
the  FIR  CO  luminosity  is  weaker  by  a  moderate  factor  than  the 
OH  and  H2  rotational  emission. 

3.2.4.  Bayesian  Analysis 

We  follow  the  Bayesian  approach  outlined  by  Ward  et  al. 
(2003)  to  quantify  the  probable  values  of  the  model  parameters. 
We  consider  an  eight-dimensional  array  of  bins  centered  on 
the  grid  points  for  which  we  have  generated  a  model  SED  and 
calculated  a  x2  value.  The  probability  that  the  actual  parameter 
set  falls  within  a  given  bin  is  proportional  to  the  product  of 
bin  size,  likelihood  L  oc  exp(— x2),  and  an  assumed  prior 
probability.  We  choose  priors  that  are  flat  in  the  logarithm  of 
each  parameter  and  that  go  to  zero  for  any  model  that  violates 
one  of  the  restrictions  listed  in  Table  2. 

In  Figure  3,  we  show  the  joint  density-temperature  proba¬ 
bility  density  functions  for  both  components,  with  contours  at 
68%,  95%,  and  99%  of  the  enclosed  probability.  The  close  sim¬ 
ilarity  between  the  95%  and  99%  (and  in  some  regions  also 
the  68%)  contours  is  due  to  a  truncation  of  the  density  func¬ 
tion  following  the  violation  of  one  of  our  model  restrictions. 
Starting  with  the  ME  component,  the  behavior  of  the  contours 
may  be  understood  as  follows.  For  either  a  fixed  Kv „  or  dv /dr, 
the  shape  of  the  CO  SED  is  approximately  conserved  if  an 
increase  in  density  is  matched  by  an  appropriate  decrease  in 
temperature.  The  pair  of  blue  curves  in  Figure  3  shows  the 
density-temperature  relation  best  fitting  the  data  for  Kvlr  =  1 
and  dv/dr  —  1000  km  s'1  pc-1.  As  the  shape  of  the  ME  SED 
is  well  constrained,  the  region  of  acceptable  parameter  space  for 
either  set  of  models  corresponds  to  a  narrow  band  centered  on 
the  best-fit  curve.  These  two  bands  intersect  at  «H2  =  109  cm'3 
and  diverge  at  lower  densities,  thereby  bracketing  a  region  of 
high  probability  filled  by  models  with  intermediate  velocity  gra¬ 
dients.  For  lower  temperatures  a  decreasing  fraction  of  the  CO 
is  excited  to  the  higher-/  states,  and  a  larger  total  mass  is  needed 
to  reproduce  the  absolute  line  fluxes.  For  «H2  >  108  cm'3  the 
gas  mass  exceeds  the  dynamical  mass,  and  consequently  a  small 
region  of  parameter  space  is  excluded.  For  higher  temperatures 
the  model  emission  does  not  fall  off  with  increasing  /  in  the 
Tupper  ~  1 8  region  as  rapidly  as  the  data  require,  and  the  lower 
quality  fits  limit  the  extent  of  the  68%  confidence  interval  above 


Figure  3.  Top:  joint  density-temperature  probability  density  function  for 
the  ME  component.  Contours  are  drawn  at  68%,  95%,  and  99%  enclosed 
probability.  The  mean  of  \og{M\{2/ Mq)  needed  to  reproduce  the  absolute  line 
fluxes  is  shown  by  the  red  curves,  and  the  logarithm  of  the  thermal  pressure 
(log(/iH2  x  Tkin)  in  K  cm-3)  is  shown  by  the  dotted  black  lines.  Blue  curves 
show  the  density-temperature  relation  giving  the  best  solution  for  Kvir  =  1  and 
dv/dr  =  103  km  s-1  pc-1  (see  the  text).  Bottom:  same  as  the  top,  but  for  the 
HE  component. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Lkin  ~  230  K.  For  l\m  >  400  and  n h2  <  10s  cm'3,  the  models 
begin  to  overproduce  the  H2  rotational  line  detections  and  upper 
limits,  and  this  accounts  for  the  sharp  cutoff  in  the  probability 
density  in  the  upper  left  region. 

The  shape  of  the  HE  SED  is  not  as  well  constrained  as 
the  ME  SED,  and  the  H2  rotational  lines  provide  a  stronger 
constraint  to  the  extent  of  the  probability  density  contours  in 
Figure  3.  The  detection  of  CO(30-29)  excludes  SEDs  peaking  at 
Lupper  <  25,  but  the  upper  limits  at  /upper  34  provide  a  softer 
restriction  on  higher  pressure  models  peaking  at  /upper  >  30. 
As  a  result,  the  dv/dr  1000  km  s'1  pc'1  restriction  no 
longer  forms  the  upper  right  boundary  of  the  density  function, 
and  models  are  allowed  to  extend  into  this  higher  density  and 
temperature  region  until  the  modeled  H2  emission  exceeds  the 
observations.  The  measured  5(1)  upper  limit  provides  the  most 
important  constraint  at  high  densities  («H2  ~  108  cm'3)  and 
low  temperatures  (Lkin  ^  300  K).  The  higher  excitation  H2 
lines  become  more  important  at  higher  temperatures,  and  the 
5(7)  transition  limits  the  upper  left  region  at  Lkin  ~  1000  K. 

In  Figure  4,  we  show  the  fully  marginalized  distribution 
functions  for  each  of  the  four  primary  model  parameters,  as  well 
as  for  A'vir  and  the  thermal  pressure  P/kg  =  nm  x  Lkjn.  The 
distribution  functions  for  the  HE  Ly n,  n m,  and  P /  kB  are  shifted 
to  higher  values  than  for  the  corresponding  ME  parameters, 
although  the  density  distributions  for  the  two  components 
contain  significant  overlap.  The  lower  limit  imposed  on  Kvn 
translates  into  a  lower  limit  to  dv/dr  that  increases  with  density 
(see  Equation  (2)).  As  such,  lower  density  solutions  are  found 
over  a  broad  range  of  velocity  gradients,  while  higher  density 
models  are  limited  to  larger  values  of  dv/dr.  This  accounts 
for  the  positive  slope  of  the  dv/dr  distribution.  Similarly,  the 
upper  limit  placed  on  dv/dr  generates  a  negative  slope  in  the 
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Figure  5.  Top:  joint  density-temperature  probability  density  function  for  the 
ME  component,  as  in  Figure  3,  but  for  various  CO  abundances.  Contours  are 
drawn  at  68%  enclosed  probability  for  [CO/H2]  =  10-5  (blue),  10-4  (black), 
and  4  x  10-4  (red).  Bottom:  same  as  the  top,  but  for  the  HE  component. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 
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Figure  4.  Probability  density  functions  for  the  four  primary  model  parameters, 
as  well  as  for  KViI  and  the  thermal  pressure  P/kg  =  hh2  x  Tjtin-  In  each  panel 
the  ME  distribution  is  shown  in  red  and  the  HE  in  blue.  The  bin  widths  are 
proportional  to  the  logarithm  of  the  model  parameter,  and  the  ordinate  indicates 
the  probability  that  the  actual  value  lies  in  the  given  bin.  The  vertical  lines 
indicate  the  median  of  each  distribution,  and  the  hatched  regions  indicate  the 
symmetric  68%  confidence  interval. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Table  3 

LVG  Model  Results 


Parameter 

ME 

HE 

Median 

68%  Range 

Median 

68%  Range 

flan  (K) 

169 

71a-347 

571 

372-896 

«H2  (cm'3) 

105  6 

104-9-106'7 

io6-4 

105'9-1071 

P/kB  (K  cm'3) 

1079 

107-4-108'7 

io9-2 

108  8— 109  8 

dv/dr  (km  s-1  pc-1) 

148 

25-1000a 

269 

71-1000a 

Ky  ir 

4.9 

la— 19 

3.9 

la-ll 

Mh2  (Mq) 

io6-7 

COs 

r 

0 

-j 

kj 

105'6 

104'9-106,3 

Note.  a  Extended  beyond  the  formal  68%  confidence  interval  to  the  truncated 
edge  of  the  distribution. 


A"vir  distribution.  As  the  distribution  functions  for  these  two 
parameters  are  heavily  influenced  by  our  prior  restrictions,  our 
modeling  produces  no  meaningful  constraint  on  the  dynamical 
state  of  the  gas. 

We  take  the  median  of  each  distribution  as  the  single  best 
estimate  of  the  parameter  value  and  indicate  this  number 
with  a  vertical  line  in  each  panel  of  Figure  4.  Most  of  the 
distribution  functions  are  relatively  symmetric,  and  we  calculate 
the  equivalent  la  uncertainties  in  the  parameter  estimation  by 


finding  the  range  of  values  symmetrically  enclosing  68%  of  the 
total  probability.  This  range  is  shown  as  the  hatched  area  under 
each  distribution.  The  results  of  this  analysis  are  summarized 
in  Table  3.  The  asymmetry  in  the  ME  7k;n  distribution  results 
from  the  limit  to  low-temperature,  high-density  parameter  space 
following  the  Kw „  is  1  and  Mgas  ML\yn  restrictions  (the  lower 
right  region  in  the  top  panel  of  Figure  3).  For  this  parameter, 
as  well  as  for  dv/dr  and  Kva,  we  extend  the  acceptable  range 
listed  in  Table  3  to  the  truncated  edge  of  the  distribution. 

3.2.5.  Varying  the  CO  Abundance 

In  the  above  analysis,  we  assumed  a  CO  abundance  of 
[CO/H2]  =  10  4,  motivated  by  abundance  measurements  in 
Galactic  molecular  clouds  (e.g.,  [CO/FF]  =  8.5  x  10~5; 
Frerking  et  al.  1982).  However,  the  CO  abundance  in  the  center 
of  NGC  1068  may  be  different.  Noting  the  general  trend  of 
higher  metallicities  in  galactic  centers,  Israel  (2009)  adopts 
an  elevated  carbon  abundance  for  the  center  of  NGC  1068 
and  derives  [CO/H2]  =  4  x  10  4.  At  the  same  time,  in 
molecular  clouds  exposed  to  intense  X-ray  fields  (as  we  expect 
for  NGC  1068;  see  Sections  5  and  6),  the  CO  abundance  may 
be  significantly  reduced  (Krolik  &  Lepp  1989;  Maloney  et  al. 
1996).  In  the  models  of  Meijerink  &  Spaans  (2005)  with  high 
ratios  of  incident  X-ray  flux  to  gas  density,  the  bulk  of  the 
gas-phase  carbon  is  not  bound  up  in  CO  until  large  depths 
(Ah  >  IO235  cm-2)  into  the  cloud.  This  yields  a  large  column 
of  warm  gas  with  reduced  CO  abundance  that  may  contribute  to 
the  FIR  CO  line  emission. 

To  explore  the  effects  of  an  altered  CO  abundance,  we  have 
repeated  our  LVG  analysis  for  a  range  of  [CO/H2]  values.  In 
Figure  5,  we  show  the  joint  density-temperature  probability 
density  functions  for  models  with  [CO/H2]  =  10  5 ,  1 0  4,  and 
4  x  1CT4.  For  a  fixed  value  of  dv/dr,  the  line  opacities  scale  as 
the  product  of  »H2  and  [CO/H2].  The  shape  of  the  CO  SED  is 
therefore  approximately  conserved  if  an  increase  in  [CO/H2] 
is  matched  by  a  decrease  in  uH2,  and  this  accounts  for  the 
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shift  to  lower  densities  for  a  larger  CO  abundance.  Increasing 
the  CO  abundance  also  leads  to  a  reduction  in  the  H2  mass 
needed  to  maintain  the  absolute  CO  fluxes  and  hence  a  smaller 
set  of  predicted  Hi  rotational  line  fluxes.  The  restriction  that 
the  modeled  Hi  emission  not  exceed  the  measured  emission 
then  provides  a  weaker  constraint  and  leads  to  an  increase  in 
the  allowed  volume  of  density-temperature  parameter  space, 
particularly  in  the  high-temperature  region.  For  [CO/FF]  = 
4  x  10”4,  the  combination  of  these  two  effects  results  in  only  a 
small  reduction  in  the  derived  density  and  H2  mass  and  a  small 
increase  in  the  derived  temperature  (Figure  5)  and  does  not 
significantly  affect  the  physical  parameter  estimates  obtained 
previously. 

If  the  CO  abundance  is  reduced  by  a  strong  X-ray  flux, 
models  indicate  that  the  hydrogen  will  also  be  largely  atomic 
(Maloney  et  al.  1996;  Meijerink  &  Spaans  2005).  We  modify  our 
LVG  calculation  to  account  for  an  atomic/molecular  mixture 
by  changing  [CO/H2]  2[CO/H],  nm  n(H)/2,  and 

Mh2  — >  M(H)  and  introducing  the  molecular  fraction  as  fmo\  — 
[H2/H].  In  each  of  these  expressions,  H  is  taken  to  represent  the 
total  number  or  mass  of  hydrogen  nuclei  in  both  atomic  and 
molecular  form.  For  collisional  excitation  of  CO  by  atomic 
hydrogen,  we  adopt  the  same  rate  coefficients  as  for  excitation 
by  H2  (see  the  discussion  in  Flower  &  Pineau  Des  Forets 
2010).  For  simplicity,  we  also  assume  that  any  reduction 
in  the  nominal  CO  abundance  is  matched  by  an  identical 
reduction  in  the  molecular  fraction  and  set  fmo\  =  2[CO/H]  x 
104.  Following  the  same  reasoning  as  discussed  above  in  the 
context  of  an  increased  CO  abundance,  a  reduction  in  the  CO 
abundance  will  increase  the  derived  density  and  total  mass. 
However,  as  the  CO-to-H2  ratio  is  conserved,  the  predicted  H2 
line  fluxes  are  to  first  order  unchanged.  The  restriction  that 
the  modeled  H2  emission  not  exceed  the  measured  emission 
then  provides  a  comparable  constraint  on  the  allowed  volume 
of  density-temperature  parameter  space  as  with  the  nominal 
CO  abundance.  For  the  ME  component,  using  2[CO/H]  = 
10~5  shifts  the  solutions  to  only  moderately  higher  densities, 
temperatures,  pressures,  and  masses  (Figure  5).  For  the  HE 
component,  the  increase  in  these  parameters  is  more  significant, 
although  the  mean  values  of  each  parameter  remain  within  the 
68%  range  for  our  nominal  [CO/H2]  =  10  4  model  (Table  3). 

4.  COMPARISON  WITH  OTHER  MOLECULAR  TRACERS 
4.1.  Physical  Parameters  and  Mass  Fraction 

The  bulk  of  the  emission  traced  by  high-resolution  millimeter 
and  near-IR  molecular  gas  maps  in  the  central  10"  arises  from 
the  ~5"  (~350  pc)  CND  (Schinnerer  et  al.  2000;  Galliano  & 
Alloin  2002;  Garcfa-Burillo  et  al.  2010).  Our  LVG  modeling 
indicates  that  no  more  than  half  of  the  CO(l-O)  emission  in  the 
CND  is  generated  by  our  ME  and  HE  components,  indicating 
that  lower  excitation  material  is  present.  The  physical  conditions 
of  this  low-excitation  component  have  been  studied  by  many 
groups.  Tacconi  et  al.  (1994)  detected  strong  CO(4-3)  emission 
toward  the  center  of  NGC  1068  and  combined  this  with  an 
interferometric  CO(l-O)  measurement  to  show  that  the  gas  was 
both  warm  (T^m  70  K)  and  dense  (nH2  )  2  x  104  cm-3). 
Subsequent  modeling  of  /Upper  C  4  transitions  of  12CO  and 
13CO  has  typically  adopted  a  fixed  Tkin  =  50  K  and  derived 
«H2  ~  104— 105  cm”3  (Sternberg  et  al.  1994;  Heifer  &  Blitz 
1995;  Usero  et  al.  2004).  Krips  et  al.  (2011)  have  obtained 
fluxes  of  the  lowest  three  transitions  of  both  12CO  and  13CO 
with  <2"  resolution.  For  a  subsection  of  the  CND  in  which 


the  gas  appears  perturbed  by  the  radio  jet,  and  hence  possibly 
shock-heated,  Krips  et  al.  (201 1 )  derive  7km  200  K  and  »H2  = 
103'5— 104-5  cm”3.  While  this  section  may  not  be  representative 
of  the  CND  as  a  whole,  this  analysis  does  suggest  that  the 
global  temperatures  may  be  higher  than  assumed  by  previous 
authors  and  similar  to  that  derived  for  our  ME  component.  We 
also  note  that  Kamenetzky  et  al.  (2011)  have  similarly  derived 
a  globally  high  temperature  ( T  >  100  K)  for  the  CND  based 
on  an  analysis  of  CS  and  other  high-density  tracers.  Given  this 
range  of  temperature  estimates,  the  clearest  difference  between 
the  ME  CO  component  and  the  lower  excitation  material  traced 
by  the  low-/  CO  lines  is  the  higher  density  («h2  ~  105'6  cm”3 
versus  «h,  <  105  cm”3).  This  suggests  a  scenario  in  which 
the  FIR  lines  trace  denser  material  in  the  CND,  which  coexists 
with  a  more  diffuse  medium  that  generates  the  millimeter  CO 
emission. 

The  mass  fraction  of  the  high  excitation  gas  may  be  estimated 
by  comparing  the  ME  and  HE  masses  derived  here  with  the 
mass  traced  by  CO(l-O).  With  a  standard  Galactic  conversion 
factor  /V(H2)//co  =  2  x  1020  cm”2  (K  km  s”1)”1,  Schinnerer 
et  al.  (2000)  estimate  a  mass  of  Mh2  =  5  x  107  M0  for  the 
CND.  This  is  similar  to  previous  estimates  from  Planesas  et  al. 
(1991,  after  correcting  their  mass  to  the  d  —  14.4  Mpc  used 
here)  and  Heifer  &  Blitz  (1995),  both  of  which  used  similar 
conversion  factors.  However,  Usero  et  al.  (2004)  have  derived  a 
much  lower  N  (H2  J//C0  =  0.3  x  1020  cm”2  (K  km  s”1)”1  (for 
the  [CO/H2]  =  10”4  used  here)  from  their  excitation  modeling 
of  the  low-/  CO  emission  from  the  CND.  Papadopoulos  & 
Seaquist  (1999b)  and  Israel  (2009)  have  also  derived  conversion 
factors  significantly  lower  than  the  Galactic  value  by  analyzing 
the  low-/  CO  emission  averaged  over  beam  sizes  of  ~L  and 
~21",  respectively.  These  authors  have  suggested  that  the  lower 
conversion  factor  arises  from  the  gas  not  being  virialized.  At  the 
same  time,  Krips  et  al.  (201 1)  have  reported  a  CO(  1-0)  flux  from 
the  central  4"  of  1 20  Jy  km  s” 1 ,  much  larger  than  the  20  Jy  km  s” 1 
value  reported  by  Schinnerer  et  al.  (2000).  For  Ico  = 

0.3  x  1020cm”2  (Kkms”1)”1  andFCo(i-<»  =  20-120  Jy  km  s”1, 
we  estimate  Mh2  =  (0.5-3)  x  107  M0.  This  is  comparable  to 
the  Mh 2  =  (0.1-5)  x  107  Mq  range  we  associate  with  our  ME 
component.  While  the  uncertainties  of  both  numbers  are  high, 
this  suggests  that  the  ME  component  makes  a  non-negligible 
contribution  to  the  total  mass  budget.  The  HE  emission  traces 
a  lower  mass  of  warmer  gas,  albeit  still  cooler  than  the  small 
amount  (M  ~  1 03  M0)  of  hot  (T  ~  2000  K)  gas  detected  in 
the  near-IR  H2  lines  (Rotaciuc  et  al.  1991;  Blietz  et  al.  1994; 
Galliano  &  Alloin  2002). 

4.2.  Line  Profiles 

4.2.1.  PACS  Lines 

Further  insight  into  the  nature  of  the  ME  and  HE  compo¬ 
nents  may  be  obtained  by  comparing  the  FIR  CO  line  profiles 
with  those  of  other  molecular  tracers.  In  particular,  we  seek  to 
understand  the  physical  origins  of  the  velocity  shift  between 
the  two  components.  To  better  demonstrate  the  spectral  shift,  in 
Figure  6,  we  show  the  composite  spectra  obtained  by  aver¬ 
aging  the  ME  (/Upper  =  14-17)  and  HE  (/upper  =  20-22, 
24,  30)  profiles.  The  PACS  spectral  resolution  changes  from 
191  to  246  km  s”1  and  127  to  311  km  s”1  over  the  wavelength 
range  corresponding  to  the  ME  and  HE  transitions,  respec¬ 
tively.  The  composite  spectra  have  been  obtained  by  smoothing 
each  line  to  311  km  s”1  resolution  (each  line  measured  with 
instrumental  resolution  Su  is  smoothed  with  a  Gaussian  kernel 
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Figure  6.  Average  profiles  of  the  ME  (red)  and  HE  (blue)  CO  lines.  Prior  to 
averaging,  each  line  is  smoothed  to  a  resolution  of  31 1  km  s-1.  The  horizontal 
error  bars  shown  underneath  the  centroid  of  each  profile  indicate  the  estimated 
drier  calibration  uncertainty  of  the  stacked  spectra. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 

of  FWHM  =  V31 12  +  Sv2)  and  resampling  to  a  common  ve¬ 
locity  grid.  Gaussian  fits  to  the  composite  profiles  yield  similar 
central  velocities  (Vme  =17  km  s_I  and  VHe  =  —64  km  s'1)  as 
obtained  from  averaging  the  central  velocities  of  the  individual 
lines  (Vme  =  17  km  s-1  and  VHe  =  —59  km  s_1;  Figure  2). 
The  uncertainties  in  the  centroids  of  the  composite  spectra  are 
dominated  by  the  spectral  calibration  uncertainties  in  the  indi¬ 
vidual  lines.  With  each  line  having  a  calibration  error  of  07 ,  we 
estimate  the  error  (er)  in  the  composite  spectra  as  a2  —  (a2) /N. 
This  gives  ctme  =  11  km  s'1  and  ohe  =  14  km  s_1,  and  we 
show  these  as  horizontal  error  bars  in  Figure  6. 

In  addition  to  CO,  the  PACS  scans  detect  strong  molecular 
emission  from  several  transitions  of  OH  and  H2O  and  weaker 
emission  from  OH+,  H20+,  and  other  molecules.  Some  of  the 
OH  and  H2O  lines  show  extended  emission  over  the  PACS  array, 
but  for  each  line  the  flux  in  the  central  spaxel  is  dominated  by 
a  compact  source  that  we  associate  with  the  CND.  The  OH+ 
and  HzO+  emission  is  also  unresolved,  although  an  association 
of  these  lines  with  the  same  molecular  gas  is  less  certain.  In 
Figure  7,  we  compare  the  central  velocities  and  widths  of  these 
lines.  For  OH,  each  point  represents  the  average  properties  of  a 
doublet,  while  for  OH+  and  H20+  each  point  is  an  average  of 
multiple  fine-structure  transitions.  On  the  right-hand  side  of  the 
top  panel,  we  show  the  mean  central  velocity  of  each  molecule, 
with  the  ME  and  HE  CO  components  separately.  The  HE  CO  is 
the  only  tracer  systematically  offset  from  the  systemic  velocity. 
In  the  bottom  panel  we  show  the  widths  of  the  ME  and  HE 
composite  profiles  as  open  points  at  X  =  109  \i  m,  where  the 
smoothed  311  km  s'  1  resolution  is  equal  to  the  instrumental 
resolution.  Overplotted  is  a  set  of  curves  showing  the  expected 
measured  line  widths  for  intrinsic  line  widths  of  0,  250,  and 
400  km  s'.  The  individual  CO  lines,  as  well  as  the  composite 
profiles,  are  consistent  with  an  intrinsic  FWHM  ~  250  km  s'1 
for  both  the  ME  and  HE  components.  The  other  molecular  lines 
are  somewhat  narrower. 

The  50-80  km  s'  1  offset  between  the  HE  lines  and  the 
other  molecular  tracers  in  Figure  7  is  <1/3  of  the  PACS 
spectral  resolution,  and  at  this  level  instrumental  effects  must 
be  considered.  The  PACS  spectral  response  depends  on  the 
illumination  of  the  slit,  such  that  a  physical  translation  of  a 
source  on  the  sky  in  the  dispersion  direction  will  produce 


rest  wavelength  [microns] 


Figure  7.  Top:  measured  line  centroids  for  the  CO,  OH  (green),  HrO  (brown), 
OH+  (purple),  and  H2CT  (cyan)  lines  detected  in  our  PACS  range  scans.  The 
CO  lines  are  color  coded  by  ME  (blue),  HE  (red),  and  intermediate  (black). 
Each  OH  point  is  an  average  of  a  doublet,  and  the  OH+  and  H20+  points 
represent  averages  over  all  detected  line-structure  lines.  The  dashed  line  shows 
the  instrumental  wavelength  shift  induced  by  pointing  offsets  of  ±2"  in  the 
dispersion  direction.  On  the  right-hand  side,  we  show  the  mean  centroid  of  each 
tracer.  Bottom:  same  as  the  top,  but  showing  the  measured  FWHM  of  each  line. 
At  A.  =  109  /tm,  we  show  the  line  widths  of  the  ME  and  HE  composite  profiles 
as  open  symbols.  Solid  curves  show  the  expected  measured  widths  for  lines 
with  intrinsic  widths  of  0,  250,  and  400  km  s  1 . 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


a  wavelength  shift  in  the  spectral  profile  (Poglitsch  et  al. 
2010).  In  the  top  panel  of  Figure  7,  we  show  the  equivalent 
velocity  shift  resulting  from  moving  a  point  source  ±2"  from 
the  center  of  the  slit.  A  comparison  with  the  CO  velocities 
shows  that  a  ~2"-3"  offset  between  the  centers  of  the  slit 
and  the  HE  CO  emission  could  be  partially  responsible  for  the 
measured  wavelength  shifts.  With  the  exception  of  CO(17-16), 
CO(24-23),  and  CO(30-29),  each  point  in  Figure  7  corresponds 
to  a  line  measured  in  the  same  set  of  concatenated  range  scans. 
Therefore,  if  all  of  the  molecular  emission  is  presumed  to  arise 
from  the  same  region  on  the  sky,  a  wavelength  shift  in  the  HE 
CO  lines  induced  by  a  pointing  offset  should  be  matched  by  a 
comparable  shift  in  the  X  1 00-150  /mi  OH  and  H2O  lines, 
and  this  is  not  observed.  Alternatively,  the  HE  CO  lines  may 
be  modeled  as  arising  from  a  region  centered  ~2"-3"  away 
from  the  source  of  the  OH  and  H2O  emission.  However,  the 
CO(30-29)  line  at  X  =  87  //  m  was  observed  separately  from  the 
series  of  concatenated  scans,  with  a  slit  position  angle  differing 
by  ~  1 80  ’.  Any  pointing-induced  shift  to  the  CO(30-29)  line 
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Figure  8.  (a)  CO(2~l)  (contours)  and  Hj  1-05(1)  (color)  images  of  the  central  4"  x  6"  from  Schinnereret  al.  (2000)  and  Muller  Sanchez  et  al.  (2009),  (b)  blue  and  red 
components  of  CO(2— 1),  with  crosses  marking  the  centers  of  the  apertures  in  panel  (a),  (c)  Hi  1-0  5(1)  image  of  the  central  0'.'4,  (e) — (i)  smoothed  Hi  1-0  5(1)  (black 
solid)  and  CO(2-l)  (black  dashed)  profiles  from  selected  apertures  in  the  Hi  ring  compared  with  the  ME  (red)  and  HE  (blue)  composite  profiles,  and  ±40  km  s-1 
horizontal  error  bar  indicating  the  calibration  uncertainty  between  the  FIR  CO  and  H2  spectra,  (d),  (j)  smoothed  Hi  1-0  S(l)  profiles  from  the  northern  and  southern 
streamer  compared  with  the  ME  and  HE  composite  profiles,  and  ±40  km  s-1  horizontal  error  bar. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


should  then  have  a  comparable  magnitude  but  be  in  the  opposite 
direction  of  the  shifts  in  the  other  HE  CO  lines,  and  this  is  also 
not  the  case.  We  therefore  argue  that  the  observed  wavelength 
shift  represents  a  real  velocity  offset,  with  instrumental  effects 
playing  no  more  than  a  minor  role. 

4.2.2.  High-resolution  Observations 

To  search  for  kinematic  substructure  in  the  CND  that  might 
explain  the  velocity  shifts  among  the  FIR  CO  lines,  we  compare 
the  FIR  line  profiles  with  those  of  high-resolution  observations 
of  CO(2-l)  (Schinnerer  et  al.  2000)  and  H2  1-0  S(  1)  (Muller 
Sanchez  et  al.  2009).  In  Figure  8(a),  we  show  maps  of  these  two 
tracers  in  the  central  4"  x  6".  At  OH  resolution,  the  CO(2-l) 


emission  from  the  CND  is  dominated  by  a  pair  of  knots  centered 
~  1"  east  and  ~  1" 5  west  of  the  AGN,  connected  by  lower  surface 
brightness  emission.  Similar  morphology  is  also  seen  in  the 
~1"  resolution  aperture  synthesis  maps  of  CN  and  SiO  (Garcla- 
Burillo  et  al.  2010)  and  of  HCN  and  HCO+  (Krips  et  al.  201 1). 
Observations  of  the  H2  1-0  S(  1 )  line  at  0'.'075  resolution  have 
resolved  the  CND  into  a  ~1"2  radius  ring  (hereafter  the  “H2 
ring”)  centered  ~0'.'6  southwest  of  the  AGN.  The  eastern  knot 
is  also  prominent  in  the  H2  map,  while  the  western  knot  is 
much  fainter.  The  H2  map  also  shows  strong  emission  from  a 
clump  ~1"  to  the  north  of  the  AGN  that  is  not  prominent  in  the 
CO(2-l)  map. 

In  panel  (b),  we  show  separate  maps  of  the  CO(2-l)  emission 
integrated  over  blue  and  red  velocities  (see  also  Figure  6  in 
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Krips  et  al.  201 1).  The  strong  emission  to  the  east  is  largely  to 
the  blue  of  the  systemic  velocity,  while  the  emission  to  the  west 
separates  into  a  blueshifted  SW  component  and  a  redshifted  NW 
component.  The  H2  velocity  field  is  generally  similar,  with  the 
exception  of  the  northern  clump  as  discussed  below.  Schinnerer 
et  al.  (2000)  model  the  CO(2-l)  kinematics  as  a  warped  disk, 
while  subsequent  study  of  the  H2  and  millimeter  tracers  has 
shown  evidence  for  additional  non-circular  motion,  possibly 
indicating  a  radially  expanding  component  (Galliano  &  Alloin 
2002;  Davies  et  al.  2008;  Garcia-Burillo  et  al.  2010;  Krips  et  al. 
201 1).  In  panels  (e)  through  (i),  we  compare  the  composite  ME 
and  HE  spectral  profiles  with  those  of  CO(2-l)  and  Hi  extracted 
from  selected  apertures.  For  each  panel,  we  have  smoothed  the 
CO(2-l)  and  H2  spectra  to  the  same  resolution  (311  km  s-1) 
and  resampled  to  the  same  grid  as  used  for  the  ME  and  HE 
composite  spectra. 

The  most  important  result  of  this  comparison  is  the  mis¬ 
match  between  the  FIR  line  profiles  and  the  broad  and  highly 
blueshifted  H2  emission  from  the  bright  knot  ~1"  to  the  north 
of  the  AGN,  shown  in  panel  (f).  In  the  high  spectral  resolu¬ 
tion  H2  map  presented  in  Galliano  &  Alloin  (2002)  this  is  the 
one  region  in  the  Hi  ring  that  displays  a  double-peaked  profile 
(the  two  peaks  become  blended  following  the  smoothing  done 
here),  with  an  extra  emission  component  to  the  blue  that  is  a 
clear  outlier  to  their  simple  rotation  plus  expansion  model.  This 
region  is  also  the  site  of  the  strongest  Br y  emission  in  the  H2 
ring,  which  displays  a  broad  (FWHM  ~  1000  km  s-1)  line.  The 
radio  jet  exits  the  nucleus  to  the  north  (Gallimore  et  al.  1996) 
and  may  interact  with  material  in  the  narrow-line  region  (Axon 
et  al.  1998).  Galliano  &  Alloin  (2002)  and  Galliano  et  al.  (2003) 
interpret  the  broad  and  complex  Bry  and  H2  line  profiles  as 
arising  from  gas  perturbed  by  the  jet  and  additionally  suggest 
that  the  ionization  resulting  from  the  jet-ISM  interaction  may  be 
responsible  for  the  strength  of  the  Bry  emission  in  this  region. 
The  smoothed  H2  profile  in  panel  (f)  is  inconsistent  with  either 
the  ME  or  HE  composite  spectrum,  and  we  conclude  that  this 
region  of  the  CND  does  not  dominate  the  PACS  CO  emission. 
In  Sections  5  and  6  we  discuss  possible  excitation  mechanisms 
for  the  FIR  CO,  including  shock  heating.  The  fact  that  we  can 
rule  out  an  origin  in  the  region  of  the  H2  ring  with  the  best 
evidence  for  molecular  gas  perturbed  by  the  jet  leads  us  to  sug¬ 
gest  that  such  jet-driven  shocks  in  the  Hi  ring  do  not  excite  the 
high-/ CO. 

Aside  from  the  bright  H2  knot  in  the  north,  the  ME  profile 
is  generally  consistent  with  much  of  the  rest  of  the  H2  ring. 
The  best  fit  is  with  the  profile  of  the  strong  H2  emission  from 
the  east.  The  profiles  from  the  NW  and  SW  are  moderately 
red-  and  blueshifted  with  respect  to  the  ME  composite,  but  a 
combination  of  these  regions,  and  indeed  of  the  H2  emission 
integrated  over  the  entire  CND  (excluding  the  northern  knot), 
would  also  generate  a  reasonable  fit.  We  note  that  the  H2  line 
widths  (FWHM  =  390-440  km  s_1)  are  larger  than  the  CO(2-l) 
line  widths  (FWHM  =  330-350  km  s-1)  in  each  panel  and 
better  match  the  widths  of  the  composite  ME  and  HE  profiles 
(FWHM  =  400-420  km  s-1).  This  may  indicate  that  the  hot  gas 
probed  by  the  H2  is  a  better  tracer  of  the  material  producing  the 
high-/  CO,  although  due  to  extinction  of  the  H2  line  (see,  e.g., 
Galliano  et  al.  2003),  the  morphology  of  the  FIR  CO  emission 
may  be  different  than  the  H2  image. 

The  blueshift  of  the  HE  emission  is  more  challenging  to 
match.  The  H2  spectra  from  the  blue  regions  in  the  E  and  S  W  are 
~65  and  ~45  km  s  1  to  the  red,  respectively,  while  the  CO(2-l) 
profiles  from  the  same  regions  are  too  narrow.  However,  the 


uncertainty  in  the  mean  centroid  of  the  HE  lines  is  ~20  km  s-1, 
while  the  H2  calibration  error  is  smaller.  Registration  errors 
between  the  velocity  frames  of  the  PACS  spectra  and  the  ground- 
based  Hi  spectra  are  also  likely  to  be  important  at  the  level  of 
~  10-20  km  s_1,  although  more  difficult  to  quantify.  With  a 
conservative  estimate  of  40  km  s'1  uncertainty  in  the  relative 
calibration  of  the  PACS  CO  and  the  H2  spectra,  the  HE  CO 
emission  only  differs  by  1.1<t-1.6<t  with  the  H2  emission  from 
the  E  or  SW  regions.  We  conclude  that  while  the  HE  emission 
profile  is  not  naturally  matched  to  the  H2  or  CO(2-l)  emission 
from  the  H2  ring,  an  association  with  the  bluest  material  to  the 
E  or  SW  may  be  within  the  measurement  errors  and  should  not 
be  excluded. 

In  panel  (c),  we  show  the  0'.'025  resolution  H2  1-0  .ST  I )  map, 
which  identifies  two  gas  clouds  streaming  toward  the  nucleus 
on  highly  elliptical  orbits  from  the  north  and  south  (hereafter 
the  “northern”  and  “southern”  streamers;  Muller  Sanchez  et  al. 
2009).  The  northern  streamer  is  connected  to  the  H2  ring  in  both 
H2  emission  and  mid-IR  continuum  and  has  been  proposed  as 
a  means  by  which  material  is  transported  to  the  AGN  (Tomono 
et  al.  2006;  Muller  Sanchez  et  al.  2009).  The  southern  streamer 
is  detected  to  within  ~10  pc  of  the  AGN  and  may  play  a  role 
in  obscuring  the  nuclear  emission.  The  southern  streamer  is 
modeled  to  lie  in  front  of  the  AGN  in  the  plane  of  the  galaxy 
and  has  a  redshifted  velocity  that  increases  from  ~50  km  s-1  to 
~80  km  s'1  as  the  gas  moves  to  within  a  projected  distance 
of  <0'.'l.  The  northern  streamer  approaches  the  AGN  from 
behind,  and  the  brightest  emission  occurs  O'.' 4  from  the  center 
at  a  blueshifted  velocity  of  30  km  s~  1 .  In  panels  (d)  and 
(j),  we  compare  the  smoothed  profiles  of  these  two  regions 
with  the  PACS  lines.  The  emission  from  the  southern  streamer 
is  detected  with  lower  S/N  but  displays  a  profile  reasonably 
consistent  with  that  of  the  ME  CO  composite.  The  line  profile 
of  the  northern  streamer  produces  an  excellent  match  to  the 
HE  composite.  The  centroids  of  the  southern  streamer  and  the 
HE  CO  composite  differ  by  ^105  km  s-1,  which  we  argue  is 
too  large  to  be  plausibly  explained  by  calibration  uncertainties 
between  the  two  data  sets.  We  conclude  that  in  addition  to  an 
origin  in  the  H2  ring  as  discussed  above,  an  origin  of  the  ME 
and  HE  components  with  the  southern  and  northern  streamers, 
respectively,  would  be  consistent  with  the  line  profiles. 

5.  HEATING  THE  ME  COMPONENT 

The  ME  CO  emission  arises  from  a  warm  (7km  ~  169  K) 
and  dense  (nHl  ~  1056  cm-3)  component,  and  with  a  total 
mass  of  Mh,  ~  10 6  7  M0  it  represents  an  important  fraction 
of  the  ISM  in  the  CND.  The  kinematic  analysis  presented  in 
Section  4.2  shows  that  this  component  may  be  readily  attributed 
to  the  H2  ring,  or  possibly  to  the  southern  streamer.  Here, 
we  consider  potential  heating  mechanisms  and  conclude  that 
X-ray-  and  shock  heating  are  both  plausible,  while  heating  by 
far-UV  photons  is  less  likely.  We  further  argue  that  no  plausible 
heating  mechanism  is  consistent  with  an  origin  in  the  southern 
streamer,  and  hence  the  emission  is  likely  associated  with  the 
H2  ring. 

5.1.  X-Ray  Heating 

The  AGN  in  NGC  1068  emits  a  hard  X-ray  luminosity  of 
Li- tokev  =  1 043 — 1 044  erg  s-1  (Iwasawa  et  al.  1997;  Colbert 
et  al.  2002).  Our  view  of  the  AGN  is  obscured  by  a  Compton- 
thick  medium,  but  the  extended  emission  detected  by  Chandra 
in  the  6-8  keV  band  demonstrates  that  the  nuclear  X-rays 
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Table  4 

Heating  Mechanisms 


ME 

HE 

Full 

XDR 

«H  =  105'75  cm-3 

Fx  =  9  erg  cm-2  s'1 
A  ~  (130  pc)2 

nn  =  105'25  cm'3 

Fx  =  160  erg  cm'2  s'1 
A  ~  (21  pc)2 

PDR 

nH  =  106,5  cm'3 

Go  =  104'75 

Lpuv  109  Lq 

nH  =  106  cm'3 
Go  =  105 
Lfuv  ~  10 10  Le 

Shock 

C-shock 

no  =  2  x  105  cm'3 
v  =  20  km  s'1 

A  ~  (150  pc)2 

C-shock 
n0  =  106  cm'3 
v  =  40  km  s'1 

A  ~  (16  pc)2 

Notes.  Details  for  the  models  used  in  Figure  9.  XDR  and  PDR  models  are  from 
Meijerink  et  al.  (2007),  ME  C-shock  model  is  from  Flower  &  Pineau  Des  Forets 
(2010),  and  HE  C-shock  model  is  from  Kaufman  &  Neufeld  (1996). 


irradiate  the  ISM  over  the  central  ~kpc  (Ogle  et  al.  2003;  Garcla- 
Burillo  et  al.  2010).  Hard  X-rays  penetrate  deeply  into  clouds 
and  efficiently  heat  large  columns  of  molecular  gas  through 
photoionization  heating  (Maloney  et  al.  1996).  The  bright  H2 
1-0  5(1)  emission  in  the  CND  has  been  attributed  to  X-ray- 
heated  gas  (Rotaciuc  et  al.  1991;  Maloney  1997;  Galliano  & 
Alloin  2002),  and  X-ray  heating  should  also  produce  strong 
emission  in  the  FIR  CO  lines  (Krolik  &  Lepp  1989). 

Meijerink  &  Spaans  (2005)  and  Meijerink  et  al.  (2007) 
present  a  detailed  photochemical  modeling  of  X-ray-dominated 
regions  (XDRs)  that  includes  predictions  for  the  emergent 
CO  line  intensities  as  a  function  of  the  gas  density  («h) 
and  incident  hard  X-ray  flux  (Fx  —  Fi- iokev)-  We  use  their 
type  A  models,  which  calculate  the  emission  from  a  parsec- 
thick  cloud  over  a  grid  covering  ;?h  =  104-106  5  cm'3  and 
Fx  —  1.6-160  erg  cm-2  s'1.  These  models  generate  CO 
line  SEDs  with  similar  shapes  as  the  isothermal  models  used 
in  our  LVG  analysis,  and  an  analogous  two-component  fit 
reproduces  the  FIR  CO  line  fluxes.  In  Figure  9(a),  we  show 
a  model  that  uses  «h  =  105375  cm'3  and  Fx  —  9  erg  cm'2  s'1 
for  the  ME  emission  (see  Table  4).  For  an  AGN  luminosity 
of  Ti-iokeV  =  1043— 1044  erg  s'1,  geometric  dilution  of  the 
radiation  field  at  the  d  ~  100  pc  distance  of  the  Hi  ring  yields 
a  flux  of  Fx  =  8.4-84  erg  cm'2  s'1,  broadly  consistent  with 
this  modeled  flux.  For  the  plane-parallel  geometry  employed 
by  Meijerink  et  al.  (2007),  the  absolute  line  luminosities  scale 
with  the  total  XDR  surface  area,  and  the  normalization  of  the 
ME  component  of  the  model  shown  in  Figure  9(a)  requires 
A  ~  (130  pc)2.  Galliano  et  al.  (2003)  model  the  Hi  ring  as  a 
section  of  a  40  pc  thick  disk.  For  a  radius  of  100  pc  the  inner 
surface  area  exposed  to  the  AGN  is  then  ~(160  pc)2,  similar 
to  the  XDR  model  requirement.  In  sum,  we  conclude  that  if  a 
substantial  fraction  of  the  H2  ring  is  exposed  to  nuclear  hard 
X-rays,  then  both  the  shape  of  the  ME  segment  of  the  CO  SED 
and  the  absolute  line  fluxes  are  naturally  reproduced  in  this  XDR 
model. 

5.2.  Far-UV Heating 

Far-UV  (FUV;  6  eV  <  hv  <  13.6  eV)  photons  offer  another 
means  of  heating  the  molecular  gas.  Photodissociation  regions 
(PDRs)  powered  by  the  FUV  radiation  from  OB  stars  in  the 
Galaxy  have  indeed  been  identified  as  prominent  sources  of 
FIR  CO  emission  (e.g.,  Kramer  et  al.  2004),  as  have  the  FUV- 
irradiated  cavities  of  protostellar  outflows  (van  Kempen  et  al. 


^upper/^B  M 


0  10  20  30  40  50 

Apper 


Figure  9.  FIR  CO  line  fluxes  and  upper  limits  measured  here,  along  with  lower-/ 
lines  from  the  literature  (see  Figure  2  for  details ),  and  XDR,  PDR,  and  shock 
model  fits,  (a)  Overplotted  is  a  two-component  XDR  fit  (solid  black),  with 
individual  components  in  red  and  blue.  Dotted  section  of  the  red  curve  shows 
an  extrapolation  of  the  model.  Upper  limits  to  the  Krolik  &  Lepp  (1989)  torus 
model  are  shown  by  the  black  (constrained  by  the  CO [44  43 1  limit)  and  green 
(constrained  by  the  stacked  limit  to  selected  /upper  =  34 — 47  transitions)  dashed 
curves,  (b)  Black  curve  shows  a  PDR  fit  to  the  full  CO  SED,  and  blue  curve 
shows  a  separate  fit  to  the  /upper  (4  1 9  lines,  (c)  Two-component  shock  fit  with 
same  color  scheme  as  panel  (a).  Blue  curve  shows  a  model  interpolated  and 
extrapolated  from  a  finite  number  of  transitions  (diamonds).  Model  parameters 
and  references  are  discussed  in  the  text  and  summarized  in  Table  4. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


2010).  However,  a  comparison  with  the  high-/  CO  emission 
from  the  prototypical  starburst  galaxy  M82  suggests  that  the 
CO  emission  from  NGC  1068  is  too  strong  to  be  powered 
by  stellar  FUV.  Herschel  observations  of  M82  have  detected 
the  submillimeter  CO  transitions  (/upper  =  4-13),  which  trace 
emission  from  gas  heated  by  shocks  (Panuzzo  et  al.  2010) 
and/or  in  PDRs  (Loenenetal.  2010).  ForanLpiR  ~  3xlOloL0 
(Telesco  &  Harper  1980;  Joy  et  al.  1987),  the  CO(13-12)/FIR 
ratio  integrated  over  the  central  43" 4  (~821  pc)  is  ~6  x  10~6. 
The  FIR  continuum  in  M82  is  produced  by  the  FUV  output 
of  young  stars,  and  we  interpret  this  CO(13-12)/FIR  ratio  as 
a  benchmark  estimate  of  the  fraction  of  FUV  radiation  that 
may  be  converted  to  FIR  CO  emission  in  an  extragalactic 
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starburst.  The  CO(14-13)/FIR  ratio  we  estimate  for  NGC  1068 
is  ~3  x  10-5,  a  factor  of  ~5  larger  than  the  CO(13-12)/FIR 
ratio  in  M82.  This  discrepancy  is  further  increased  by  noting  that 
only  a  minor  fraction  of  the  FIR  continuum  from  the  CND  in 
NGC  1068  originates  in  young  stars.  The  nuclear  stellar  cluster 
in  NGC  1068  has  a  characteristic  age  of  200-300  Myr,  and  hence 
Choi/Tfc  70  (for  a  starburst  with  exponential  decay  timescale 
of  tsf  =  100  Myr;  Davies  et  al.  2007).  For  the  LK  measured 
by  Davies  et  al.  (2007),  this  corresponds  to  Lhni  <  3  x  109  Lq, 
more  than  five  times  lower  than  the  Lfir  measured  here.  PDRs 
heated  by  the  output  of  this  cluster  would  therefore  have  to 
be  ~25  times  more  efficient  in  converting  the  incident  FUV 
radiation  to  FIR  CO  emission  than  the  PDRs  in  the  nucleus 
of  M82.  Loenen  et  al.  (2010)  model  the  /upper  =  1-13  CO 
emission  from  M82  with  a  three-component  PDR  model,  in 
which  the  CO(13-12)  emission  is  almost  entirely  generated  by 
the  highest  excitation  component.  The  CO(13-12)/FIR  ratio  of 
this  component  is  comparable  to  the  measured  CO(14-13)/FIR 
ratio  in  NGC  1068,  but  the  total  modeled  ratio  is  a  factor  of  ~  14 
lower  due  to  the  FIR  emission  from  the  lower  excitation  material. 
Increasing  the  net  ratio  by  a  factor  of  ~25  would  require  that 
the  total  FIR  be  dominated  by  the  highest  excitation  component 
and  hence  a  significant  suppression  of  the  lower  excitation  PDR 
emission.  We  consider  such  an  extreme  repartitioning  of  the 
stellar  FUV  between  M82  and  NGC  1068  less  likely  than  the 
XDR  scenario  discussed  in  Section  5.1. 

The  CO  emission  is  more  plausibly  attributed  to  PDRs 
powered  by  the  AGN  FUV  radiation,  and  in  many  ways  this  is 
an  attractive  option.  The  intrinsic  AGN  luminosity  in  NGC  1068 
of  Tfuv  ^  7.4  x  1043  erg  s  1  (Pier  et  al.  1994)  is  similar  to 
the  observed  Lfir,  supporting  the  idea  that  the  FIR  continuum 
is  emission  from  dust  grains  heated  by  the  AGN.  The  measured 
60  jtm/100  jim  flux  ratio  of  ~1.3  implies  that  the  incident 
FUV  flux  on  these  grains  corresponds  to  Go  ~  105  (where 
Lfuv  =  Go  x  1.6  x  10~3  erg  s-1  cm'2)  (Abel  et  al.  2009), 
which  would  also  be  expected  if  the  AGN  LFuv  is  absorbed 
at  the  d  ~  100  pc  distance  of  the  FU  ring.  To  estimate  the 
CO  emission  produced  in  these  PDRs,  we  employ  the  type  A 
PDR  models  of  Meijerink  et  al.  (2007),  which  adopt  the  same 
geometry  and  gas  densities  as  the  XDR  models  and  span  an  FUV 
intensity  range  of  Go  =  102-105.  In  dense  PDRs  most  of  the  CO 
forms  at  cloud  depths  greater  than  Ay  &  2  at  gas  temperatures 
T  <  100  K,  but  a  secondary  CO  abundance  peak  at  Ay  ~  0.6 
with  T  ~  800  K  follows  from  a  local  enhancement  of  OFI 
(Sternberg  &  Dalgarno  1995).  The  CO  line  SEDs  generated 
in  the  Meijerink  et  al.  (2007)  models  are  characterized  by  two 
distinct  peaks,  which  we  associate  with  these  warm  and  hot  CO 
phases.  In  Figure  9(b),  we  show  a  model  with  nH  =  106  cm'2 
and  Go  =  10s  that  produces  an  intriguing  match  to  most  (with 
the  exception  of  CO[30-29])  of  the  FIR  CO  line  SED  (solid 
black  line),  with  the  warm  and  hot  phases  reproducing  the  ME 
and  HE  components,  respectively.  Assuming  that  the  emitted 
FIR  continuum  is  equal  to  the  incident  FUV,  this  model  predicts 
a  CO(14-13)/FIR  ratio  of  4  x  10-5.  This  is  similar  to  the 
observed  ratio  in  high-Go  Galactic  PDRs  (Kramer  et  al.  2004) 
and  the  value  measured  here  for  NGC  1068. 

However,  there  are  at  least  two  reasons  to  reject  this  scenario. 
First,  the  AGN  UV/optical  emission  escapes  to  d  ~  100  pc 
distances  only  in  the  ionization  cone  (P.A.  15°;  Macchetto 

et  al.  1994),  while  the  brightest  emission  from  the  Hi  ring  is 
to  the  east.  If  the  warm  molecular  gas  traced  by  the  ME  CO 
emission  is  assumed  to  be  cospatial  with  the  hot  gas  traced  by 
the  H2  1-0  S(  1)  line,  then  we  would  expect  little  interaction 


between  this  warm  gas  and  the  nuclear  FUV.  Additionally,  any 
single  PDR  model  that  simultaneously  matches  the  ME  and  HE 
emission  is  inconsistent  with  the  kinematic  evidence  that  these 
two  sets  of  lines  are  tracing  physically  distinct  components.  In 
general,  we  require  that  any  model  reproducing  the  Lupper  <17 
lines  must  underpredict  the  fluxes  in  the  JuppL,r  >  20  transitions. 
In  contrast,  all  of  the  models  produced  by  Meijerink  et  al.  (2007) 
that  fit  the  ME  section  of  the  CO  SED  either  match  or  exceed 
the  observed  HE  line  fluxes.  Consequently,  we  exclude  PDRs  as 
a  potential  source  of  the  ME  emission.  We  stress  that  the  latter 
argument  depends  sensitively  on  an  accurate  modeling  of  the 
7upPei  >  20  CO  emission  from  the  hot  surfaces  of  PDRs.  Future 
Herschel- PACS  studies  of  the  FIR  CO  line  SEDs  in  Galactic 
PDR  templates  will  be  useful  in  further  evaluating  this  PDR 
scenario. 


5.3.  Shock  Heating 

Shock  heating  offers  a  simple  means  of  exciting  the  high- 
J  CO  emission  and  is  generally  the  preferred  mechanism  for 
producing  the  highest  excitation  molecular  gas  in  Galactic 
sources  (e.g.,  Sempere  et  al.  2000).  Using  the  C-shock  models 
of  Flower  &  Pineau  Des  Forets  (2010),  we  find  that  the  ME 
component  of  the  CO  SED  can  be  fit  with  a  pre-shock  density  of 
«H  =  2  x  10s  cm'3  and  shock  velocities  of  vs  —  10-20  km  s_1 . 
The  contributions  of  H2O  and  H2  to  the  total  cooling  budget 
increase  rapidly  with  shock  velocity  in  these  models,  and 
keeping  vs  20  km  s' 1  is  necessary  to  prevent  the  predicted 
H2O  and  H2  line  fluxes  from  exceeding  the  measured  values. 
Decreasing  the  shock  velocity  from  vs  =  20  to  10  km  s'  1 
generates  weaker  CO  lines  and  requires  increasing  the  total  cross 
section  from  A  ~  (150  pc)2  to  (210  pc)2  to  match  the  absolute 
line  fluxes.  In  Figure  9(c),  we  combine  the  vs  =  20  km  s'  1 
model  with  a  separate  C-shock  fit  to  the  HE  CO  emission. 

What  shock  mechanism  could  produce  the  ME  CO  line 
emission?  Mechanical  heating  from  stellar  feedback  has  been 
proposed  as  the  energy  source  behind  the  submillimeter  CO 
emission  in  the  M82  starburst  (Panuzzo  et  al.  2010)  and  in 
other  galactic  nuclei  (Hailey-Dunsheath  et  al.  2008;  Nikola  et  al. 
2011),  but  the  contrast  in  CO/FIR  ratios  between  NGC  1068 
and  M82  discussed  in  Section  5.2  argues  against  a  similar 
mechanism  here.  Furthermore,  we  can  use  the  results  of  Davies 
et  al.  (2007)  to  estimate  the  mechanical  power  injected  into  the 
ISM  by  SNe  in  the  nuclear  cluster  of  NGC  1068.  These  authors 
model  the  SN  rate  (SNR)  to  LK  ratio  as  <6  x  10-11  yr' 1  L/j1  in 
this  cluster,  which  combined  with  their  measured  Lk  yields 
SNR  <  2.4  x  1 0''.  Following  Loenen  et  al.  (2008)  and 
estimating  a  mechanical  energy  release  of  1051  erg  per  SN, 
with  10%  dissipated  in  molecular  gas,  yields  a  heating  rate  of 
<2  x  106  Lq.  This  is  at  least  a  factor  of  ~9-37  times  lower 
than  the  mechanical  luminosity  L  =  1/2 pu3A  required  by 
the  shock  models  discussed  above.  Jet-ISM  interactions  are 
another  potential  source  of  shock  heating  in  Seyfert  nuclei,  but 
as  discussed  in  Section  4.2,  the  mismatch  in  line  profiles  between 
the  FIR  CO  lines  and  the  disturbed  H2  1-0  S{  1 )  profile  ~  1"  north 
of  the  AGN  suggests  that  jet-driven  shocks  in  the  H2  ring  may 
not  be  important. 

Alternatively,  we  note  that  the  cross  sections  required  to 
normalize  these  plane-parallel  shock  models  are  similar  to  the 
estimated  cross  section  of  the  H2  ring  as  viewed  from  the  center. 
The  kinematics  drawn  from  the  high-resolution  H2  1-0  5(1)  and 
millimeter- wave  molecular  gas  maps  require  a  radial  expansion 
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component  to  the  H2  ring  (Galliano  &  Alloin  2002;  Davies 
et  al.  2008;  Krips  et  al.  201 1).  Galliano  &  Alloin  (2002)  model 
the  H2  dynamics  by  combining  a  rotational  component  with  a 
v  =  140  km  s~  1  radially  expanding  component  that  generates 
1  /3  of  the  total  line  emission,  and  Krips  et  al.  (201 1)  construct 
a  similar  model  to  explain  the  CO  dynamics.  We  suggest 
that  the  interaction  of  the  outflowing  gas  with  non-outflowing 
material  in  the  H2  ring  offers  the  most  plausible  source  of 
shock  heating.  We  recall  that  our  excitation  model  requires  the 
dense  gas  associated  with  the  ME  CO  emission  to  be  mixed 
with  lower  density  (n H2  105  cm~3)  material  responsible  for 
the  millimeter-wave  CO  emission  (Section  4.1).  Assuming  that 
the  nv2  product  is  conserved  for  shocks  propagating  through 
an  inhomogeneous  medium  (Klein  et  al.  1994),  the  moderate 
( v  <  20  km  s_1)  velocities  we  require  could  ultimately  be 
produced  in  i>  ~  140  km  s”1  shocks  triggered  in  lower  density 
gas. 

5.4.  Southern  Streamer 

In  Section  4.2,  we  showed  that  the  ME  CO  line  profile  may 
be  consistent  with  that  of  the  H2  1-0  5(1)  emission  from  the 
southern  streamer.  However,  such  an  association  is  inconsistent 
with  either  the  X-ray-  or  shock-heating  scenarios  outlined 
above.  In  both  of  these  models  the  required  cross  section  of 
A  ~  (130  pc)2-(150  pc)2  is  significantly  larger  than  the  ~20  pc 
size  of  the  southern  streamer.  Equivalently,  the  absolute  intensity 
of  the  CO  emission  from  this  cloud  would  have  to  be  more  than 
an  order  of  magnitude  larger  than  predicted  by  the  models.  As 
such,  we  argue  that  the  ME  CO  emission  does  not  arise  from 
the  southern  streamer,  but  most  likely  from  the  H2  ring. 

5.5.  Summary  and  Discussion 

In  summary,  we  conclude  that  the  ME  CO  emission  arises 
from  either  X-ray-  or  shock-heated  gas  in  the  H2  ring.  The 
challenge  of  unambiguously  determining  the  heat  source  for  this 
gas  is  similar  to  the  situation  for  the  warm  and  hot  gas  traced  by 
the  H2  rotational  and  rovibrational  lines  in  NGC  1068  (Rotaciuc 
et  al.  1991;  Lutz  et  al.  2000)  and  in  larger  samples  of  Seyferts 
(Davies  et  al.  2005;  Rodrfguez-Ardila  et  al.  2005;  Roussel  et  al. 
2007).  The  link  between  the  ME  CO  and  the  H2  emission  in 
NGC  1068  is  not  immediately  clear,  although  at  least  to  within 
the  PACS  resolution  the  similarity  of  the  ME  CO  and  H2  1-0  5(1) 
line  profiles  suggests  a  connection  (Figure  8).  Additionally,  we 
note  that  the  H2  1-0  5(1)  brightness  is  quantitatively  reproduced 
with  a  similar  XDR  model  as  used  here  for  the  ME  CO  (although 
using  a  lower  density  of  n  =  105  cm  3 ;  Maloney  1 997 ;  Galliano 
et  al.  2003),  while  the  shock  model  shown  in  Figure  9(c)  also 
produces  an  H2  1-0  5(1)  flux  within  a  factor  of  ~2  of  the  total 
CND  emission.  Further  joint  modeling  of  the  CO,  H2,  and  other 
molecular  emission  in  NGC  1068,  as  well  as  FIR  CO  data  on 
a  larger  sample  of  comparison  sources,  will  be  useful  in  better 
understanding  the  nature  of  the  ME  CO  component. 

Here  we  offer  two  reasons  for  preferring  the  X-ray-heating 
scenario.  First,  the  hard  X-ray  luminosity  of  NGC  1068  is 
reasonably  well  established  either  through  an  analysis  of  the 
directly  observed  (scattered)  emission  or  by  applying  scaling 
relations  established  for  type  1  systems  (Iwasawa  et  al.  1997; 
Colbert  et  al.  2002).  As  discussed  above,  combining  this 
luminosity  with  reasonable  estimates  of  the  gas  density  and 
the  CND  geometry  naturally  produces  the  ME  CO  line  SED.  In 
contrast,  it  is  not  clear  whether  the  shocks  in  the  CND  are 
sufficient  to  dissipate  enough  mechanical  energy  at  the  low 


velocities  needed  to  power  the  CO.  As  we  argued  above,  jet- 
driven  shocks  in  the  H2  ring  do  not  provide  a  good  fit,  while 
the  energetics  of  the  possible  shocks  arising  from  the  radial 
expansion  of  the  H2  ring  have  yet  to  be  demonstrated.  Second, 
a  number  of  authors  have  noted  that  the  chemical  composition 
of  the  CND  is  best  described  through  a  model  of  X-ray-driven 
chemistry  (Usero  et  al.  2004;  Krips  et  al.  2008,  2011;  Garcfa- 
Burillo  et  al.  2010).  The  OH+  and  H20+  emission  detected  in 
our  PACS  scans  offers  further  evidence  (van  der  Werf  et  al. 
2010;  Rangwala  et  al.  2011)  and  will  be  discussed  in  a  future 
paper.  If  the  nuclear  X-rays  are  responsible  for  the  anomalous 
molecular  abundances  in  the  CND,  it  is  likely  that  they  also  play 
an  important  role  in  the  energetics. 

6.  HEATING  THE  HE  COMPONENT 

The  HE  CO  emission  arises  from  a  small  mass  (Mh2  ~ 
105-6  M0)ofwarm(7kin  ~  571  K)  and  dense  (nH2  ~  106  4cm_3) 
material  that  represents  only  a  minor  fraction  (<  1 0%)  of  the  total 
gas  in  the  CND.  The  kinematic  analysis  presented  in  Section  4.2 
suggests  that  this  component  may  potentially  be  associated  with 
the  most  blueshifted  emission  in  the  east  or  west  of  the  H2  ring, 
or  with  the  clump  of  infalling  gas  ~0'.'4  north  of  the  AGN. 
Here,  we  consider  potential  heating  mechanisms  and  argue  that 
the  HE  CO  emission  arises  from  either  X-ray-heated  gas  in  the 
northern  streamer  or  shock-heated  material  in  either  the  northern 
streamer  or  H2  ring. 

6.1.  X-Ray  Heating 

In  Section  5.1,  we  used  the  Meijerink  et  al.  (2007)  XDR 
models  to  generate  a  two-component  fit  to  the  FIR  CO  emission 
and  overplotted  a  sample  solution  in  Figure  9(a).  For  the 
mh  =  105—106  5  cm~3  density  range  the  ME  component  requires 
Fx  =  5.1-16  erg  cm-2  s  ,  while  for  the  same  densities 
the  HE  component  requires  the  maximum  (or  higher)  Fx  = 
160  erg  cm~2  s~  1  used  in  the  Meijerink  et  al.  (2007)  model 
grid.  An  important  outcome  of  this  XDR  modeling  is  therefore 
that  the  HE  component  requires  irradiation  by  an  X-ray  field  at 
least  an  order  of  magnitude  stronger  than  the  ME  component. 
This  is  difficult  to  achieve  if  both  components  are  presumed  to 
arise  from  the  H2  ring.  The  eastern  segment  of  the  ring  lies  no 
more  than  a  factor  of  ~  1.5  closer  to  the  AGN  than  the  western 
segment,  so  a  variation  in  the  distance  to  the  AGN  across  the 
ring  appears  unlikely  to  produce  a  factor  of  ~10  variation  in 
the  radiation  field  strength.  Garcfa-Burillo  et  al.  (2010)  have 
used  the  6-8  keV  emission  detected  with  Chandra  to  trace 
the  penetration  of  nuclear  hard  X-rays  into  the  cold  neutral 
ISM.  Their  image  suggests  a  relatively  isotropic  illumination 
of  the  central  few  hundred  parsecs  (see  also  Ogle  et  al.  2003). 
Specifically,  while  the  hard  X-rays  may  be  expected  to  escape  to 
~100  pc  scales  more  easily  in  the  ionization  cone,  the  6-8  keV 
band  image  is  not  significantly  enhanced  along  this  direction 
within  the  H2  ring. 

The  most  straightforward  way  of  generating  a  higher  X-ray 
intensity  is  to  attribute  the  HE  emission  to  the  northern  streamer. 
In  projection,  the  northern  streamer  lies  a  factor  of  ~3  closer 
to  the  AGN  than  does  the  H2  ring  (~0'.'4  versus  ~  172;  Muller 
Sanchez  et  al.  2009)  and  is  therefore  expected  to  see  a  ~9  times 
stronger  Fx.  A  fit  to  the  HE  component  with  nH  =  105  25  cm  3 
and  Fx  —  160  erg  cm~2  s_1  (Figure  9(a))  requires  a  surface  area 
of  ~(21  pc)2,  close  to  the  ~14  pc  (~0'.'2)  projected  lateral  size 
of  this  clump.  We  therefore  argue  that  if  X-ray  heating  accounts 
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Figure  10.  Ratio  of  CO(19-18)  emission  from  a  PDR  to  an  XDR  (black),  using 
the  models  of  Meijerink  et  al.  (2007).  The  incident  X-ray  flux  in  the  XDR  model 
is  six  times  weaker  than  the  FUV  flux  in  the  PDR  model.  The  fraction  of  the 
measured  Lfir  that  must  be  attributed  to  this  PDR  is  shown  in  red.  The  region 
of  high  density  and  high  Go  that  can  match  the  HE  CO  line  SED  (excluding 
CO[30-29])  without  overpredicting  the  lower-7  fluxes  is  indicated  with  the 
dotted  lines. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 

for  both  the  ME  and  HE  components,  the  HE  component  is  most 
plausibly  associated  with  this  infalling  gas. 

6.2.  Far-UV Heating 

Using  the  Meijerink  et  al.  (2007)  PDR  models,  we  find  that 
for  a  proper  choice  of  high  density  (nH  >  106  cm  3)  and  strong 
FUV  field  (Go  >  104),  the  models  reproduce  the  7upper  ~  19-24 
lines  while  underpredicting  the  lower-7  fluxes.  An  example 
model  with  «h  =  106'5  cm-3  and  Go  =  104'75  is  shown  in 
Figure  9(b)  (thin  blue  line).  The  FUV  continuum  associated 
with  this  component  is  Efuv  ~  2  x  109  L0,  a  factor  of  ~10 
less  than  the  observed  Lfir-  Thus,  we  may  consider  a  scenario 
in  which  the  ME  CO  emission  arises  from  X-ray-  or  shock- 
heated  gas  in  the  H2  ring,  while  the  bulk  of  the  HE  emission 
originates  in  a  trace  amount  of  PDR  material  that  makes  only 
a  minor  contribution  to  the  total  continuum  emission.  These 
PDR  models  all  fail  to  reproduce  the  CO(30-29)  transition, 
however,  which  in  this  scenario  must  therefore  arise  from  a 
third  component  of  yet  more  highly  excited  gas. 

Any  gas  in  the  CND  exposed  to  the  FUV  from  the  AGN  must 
also  be  irradiated  by  hard  X-rays.  For  an  FUV  origin  of  the  HE 
CO  lines,  the  physical  properties  and  radiative  environment  of 
the  clouds  must  therefore  be  such  that  the  FUV  is  more  effective 
than  the  X-rays  in  generating  7UpPer  >19  CO  emission.  Under 
what  conditions  would  this  occur?  A  first  approach  to  addressing 
this  question  is  to  separately  consider  the  CO  emission  from  a 
PDR  with  that  of  a  properly  selected  XDR.  In  Figure  10,  we 
compare  the  CO(  19-18)  (the  highest-7  transition  calculated  over 
the  full  model  grid)  flux  from  the  Meijerink  et  al.  (2007)  PDR 
models  with  that  of  a  Meijerink  et  al.  (2007)  model  XDR  with 
s»l/6  of  the  incident  flux.  This  flux  ratio  is  chosen  to  match 
the  intrinsic  LFUv/ix  ~  6  ratio  of  the  AGN  in  NGC  1068 
(Pier  et  al.  1994)  and  is  therefore  appropriate  for  a  cloud  seeing 
the  unattenuated  AGN  emission.  The  red  curves  in  Figure  10 
show  the  fraction  of  the  observed  FIR  continuum  (assuming 
LFUV  =  LFIR)  that  must  arise  from  the  model  PDR  in  order  to 
account  for  the  absolute  CO(20-19)  line  flux.  The  region  in  the 
upper  right  quadrant  enclosed  by  the  dotted  lines  indicates  the 


subset  of  high  excitation  PDR  models  that  reproduce  the  HE  CO 
lines  (excluding  CO[30-29])  while  underpredicting  the  lower-7 
fluxes — a  requirement  for  this  solution  given  the  velocity  shifts 
between  the  ME  and  HE  lines.  For  clouds  with  «h  >  106  cm”3, 
Go  >  104,  and  Go/wh  /5  10”L5  cm3,  this  approach  suggests 
that  the  FUV  is  more  important  than  the  X-rays  in  generating 
Supper  >  19  CO  emission,  while  only  a  minor  fraction  of  the 
AGN  FUV  luminosity  would  be  required  to  power  such  a  PDR. 

Where  in  the  nuclear  region  might  these  conditions  be 
satisfied?  The  UV/optical  emission  from  the  AGN  escapes  to 
large  distances  within  the  ionization  cone,  which  runs  roughly 
north-south  in  projection.  If  the  HE  CO  emission  arises  from 
the  H2  ring,  however,  the  line  blueshifts  would  suggest  an  origin 
to  the  east  or  west  (Section  4.2),  where  little  nuclear  FUV 
penetrates.  The  northern  streamer  may  have  a  direct  view  of 
the  AGN,  which  at  a  distance  of  d  &  40  pc  would  correspond 
to  Go  ~  105'4.  However,  the  covering  factor  of  this  material 
(for  a  size  of  ~14  pc)  is  only  ~0.01.  The  contents  of  Figure  10 
indicate  that  achieving  the  absolute  HE  CO  line  fluxes  with 
such  a  small  covering  factor  would  require  densities  far  larger 
than  the  «H2  ~  106'4  cm”3  indicated  by  our  LVG  modeling 
(Section  3.2).  We  therefore  argue  that  the  HE  CO  emission  is 
unlikely  to  be  FUV-powered,  although  more  detailed  modeling 
of  the  conversion  of  FUV  radiation  to  high-7  CO  at  high  densities 
(«h  >  106'5  cm”3)  and  FUV  intensities  (Go  >  10s),  and  in  the 
presence  of  an  additional  X-ray  field,  would  be  useful  for  a  more 
quantitative  evaluation  of  this  scenario. 

6.3.  Shock  Heating 

The  HE  CO  component  is  fit  by  a  C-shock  model  with  pre¬ 
shock  density  «H2  —  106  cm”3,  velocity  vs  —  40  km  s”1, 
and  cross  section  A  ~  (16  pc)2  (Kaufman  &  Neufeld  1996) 
(Figure  9(c)).  As  discussed  in  Sections  4.2  and  5.3,  jet-ISM 
interactions  in  the  H2  ring  are  unlikely  to  produce  the  FIR  CO 
emission,  but  a  jet-induced  shock  in  the  northern  streamer  is 
a  more  plausible  source  of  the  HE  CO  lines.  The  radio  jet 
changes  direction  in  the  vicinity  of  the  brightest  clump  in  the 
northern  streamer,  evidence  that  the  fast-moving,  low-density 
material  in  the  jet  is  colliding  with  and  being  diverted  by 
the  dense  molecular  gas  (Muller  Sanchez  et  al.  2009).  This 
interaction  should  produce  shocks  in  the  molecular  material,  and 
the  presence  of  H2O  masers  in  this  cloud  supports  this  picture 
(Gallimore  et  al.  2004).  The  modeled  A  ~  (16  pc)2  cross  section 
matches  the  ~14  pc  size  of  this  clump,  consistent  with  this 
scenario.  The  jet  mechanical  power  estimated  in  the  bow  shock 
model  of  Wilson  &  Ulvestad  (1987)  is  ~2  x  108  L0,  a  factor 
of  ~3  larger  than  the  mechanical  luminosity  L  =  1/2 pv^A 
of  the  shock  model  discussed  here.  Jet-induced  shocks  are 
therefore  energetically  feasible,  although  they  would  require 
an  efficient  mechanism  for  converting  kinetic  energy  at  high 
velocities  into  slow  molecular  shocks.  Alternatively,  shocks  in 
the  bluest  regions  of  the  H2  ring,  possibly  associated  with  the 
ring  expansion,  could  also  generate  the  HE  CO. 

An  important  constraint  on  this  shock  modeling  is  the 
need  to  match  the  CO  brightness  while  not  overproducing  the 
rovibrational  H2  emission.  For  the  C-shock  model  discussed 
above  the  predicted  H2  1-0  5(1)  flux  is  a  factor  of  ~  1 90  larger 
than  observed  in  the  northern  streamer  (Muller  Sanchez  et  al. 
2009)  and  ~19  larger  than  in  the  bright  eastern  clump  in  the  H2 
ring  (Galliano  &  Alloin  2002).  In  this  region  of  the  Kaufman 
&  Neufeld  (1996)  model  parameter  space,  the  H2/CO  intensity 
ratio  is  a  steep  function  of  shock  velocity,  while  the  CO  emission 
is  sensitive  to  both  the  velocity  and  pre-shock  density.  Finely 
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tuning  vs  and  n o  (going  to  lower  vs  and  higher  hq)  may  therefore 
yield  a  solution  that  reduces  the  H2  emission  while  conserving 
the  CO  line  SED.  Extinction  of  the  2  /rm  EE  emission  may 
also  be  important,  particularly  in  the  northern  streamer.  Muller 
Sanchez  et  al.  (2009)  estimate  a  column  of  Ah  ~  8  x  1024  cm”2 
in  the  southern  streamer.  Assuming  a  similar  value  for  the 
northern  streamer,  and  for  Ak/Nh  =  (2  x  1022  cm”2)”1, 
the  resulting  Ak  ~  400  is  more  than  sufficient  to  provide 
the  necessary  attenuation  in  a  mixed  dust  model.  A  /-shock 
model  for  the  HE  CO  would  also  generate  much  weaker  H2 
emission  (Hollenbach  &  McKee  1989),  although  the  required 
cross  section  would  then  be  a  factor  of  ~9  higher  than  for  a 
C-shock  and  no  longer  match  the  northern  streamer  size.  In 
sum,  a  number  of  scenarios  may  be  envisioned  to  attribute  the 
HE  CO  lines  to  a  shock  in  the  H2  ring  or  the  northern  streamer, 
while  also  satisfying  the  H2  1-0  S(l)  brightness. 

6.4.  Summary  and  Discussion 

We  conclude  that  the  HE  component  most  plausibly  arises 
from  X-ray-  or  shock-heated  gas  in  the  northern  streamer  or 
shock-heated  gas  in  the  H2  ring.  For  the  ME  component,  we 
argued  in  Section  5.5  that  energetic  considerations  and  the 
ancillary  evidence  for  an  XDR  chemistry  favored  X-ray  over 
shock  heating,  but  the  situation  for  the  HE  component  is  less 
clear.  The  amount  of  mechanical  power  available  in  the  jet 
or  other  sources  is  indeed  uncertain,  but  the  strong  evidence 
for  a  jet  interaction  with  the  northern  streamer  suggests  that 
jet-induced  shocks  in  this  cloud  be  given  full  consideration. 
Additionally,  the  OH+  and  H20+  emission  lines  that  pinpoint 
an  XDR  chemistry  are  emitted  close  to  the  galaxy  systemic 
velocity,  with  no  detected  emission  at  the  blueshifted  velocity 
of  the  HE  CO  lines  (Figure  7).  Thus,  while  nuclear  X-rays  may 
dominate  the  energetics  and  chemistry  in  the  CND,  an  origin  of 
the  HE  CO  lines  in  a  small  amount  of  shock-heated  gas  must 
also  be  considered. 

7.  CONSTRAINTS  ON  THE  NUCLEAR  OBSCURATION 

7.1.  Could  the  Detected  Emission  Arise  from 
Compton-thick  Gas  near  the  AGN? 

The  hard  X-rays  emitted  by  the  AGN  in  NGC  1068  are 
obscured  by  a  Compton-thick  medium,  which  may  have  a 
column  density  as  high  as  /Vh  >  1025  cm”2  (Matt  et  al.  1997). 
Such  a  high  column  should  provide  enough  X-ray  shielding 
to  enable  the  formation  of  CO  and  other  molecules,  and  this 
gas  may  be  sufficiently  warm  and  dense  to  excite  the  FIR 
CO  transitions  (Krolik  &  Lepp  1989).  Molecular  observations 
of  NGC  1068  have  identified  two  possible  candidates  for 
this  obscuring  material  within  the  central  <20  pc.  Radio 
observations  at  22  GHz  have  detected  a  string  of  H2O  masers 
that  appear  to  trace  a  thin,  rotating  disk,  centered  on  the  AGN, 
with  inner  and  outer  radii  of  ~0.65  and  ~1.1  pc,  respectively 
(Gallimore  et  al.  2004,  and  references  therein).  At  larger 
distances,  the  southern  streamer  is  observed  to  approach  to 
within  ~10  pc  of  the  AGN  and  may  be  the  outer  part  of  an 
amorphous,  clumpy  structure  obscuring  the  nucleus  (Muller 
Sanchez  et  al.  2009).  Throughout  this  paper,  we  have  attributed 
the  detected  FIR  CO  emission  to  gas  associated  with  either  the 
H2  ring  or  the  northern  streamer,  with  the  implication  that  we 
are  not  detecting  the  obscuring  medium.  In  Sections  4.2.2  and 
5.4,  we  explicitly  argued  that  the  line  profiles  and/or  absolute 
intensities  of  the  ME  and  HE  components  are  inconsistent  with 


either  component  arising  from  the  southern  streamer.  However, 
is  it  possible  that  the  HE  emission  arises  from  the  maser  disk, 
or  elsewhere  within  the  central  few  parsecs?  And  could  this  gas 
provide  the  Ah  ~  1025  cm-2  obscuring  column? 

The  geometrical  constraints  imposed  by  the  LVG  modeling 
do,  in  fact,  allow  for  the  HE  component  to  comprise  a  parsec- 
scale,  Compton-thick  structure.  The  LVG  model  used  here 
assumes  emission  from  a  collection  of  spherical  clouds  and 
retains  the  freedom  to  arrange  these  clouds  in  an  arbitrary 
configuration.  As  an  example,  we  consider  smoothly  distributing 
the  gas  in  a  spherical  shell  with  a  finite  covering  factor  (/cov). 
Hard  X-ray  surveys  conducted  with  International  Gamma- 
Ray  Astrophysics  Laboratory /IBIS  (Malizia  et  al.  2009)  and 
Swift- Burst  Alert  Telescope  (Burlon  et  al.  2011)  indicate  that 
~20%-25%  of  AGNs  are  Compton-thick,  and  in  the  spirit  of 
unification  we  therefore  adopt  /cov  =  0.25.  We  require  the  shell 
to  have  a  column  density  of  Ah  =  1025  cm'2,  in  which  case 
a  choice  of  density  and  total  mass  determines  the  inner  radius 
(Rm).  For  the  set  of  good-fitting  LVG  solutions  discussed  in 
Section  3.2.3,  we  find  values  of  Rm  =  0-7  pc.  This  range  of 
radii  allows  structures  matched  in  size  to  the  maser  disk,  as  well 
as  moderately  more  extended  configurations. 

The  ~250  km  s'1  FWHM  of  the  FIR  CO  lines  is  considerably 
lower  than  the  ~600  km  s'1  velocity  range  of  the  maser 
spots.  If  the  HE  component  is  attributed  to  a  uniform  disk, 
virial  considerations  then  require  a  larger  size  scale  than  the 
r  =  0.65-1.1  pc  radial  extent  of  the  maser  disk.  Translating  the 
observed  line  width  to  a  disk  size  is  beyond  the  scope  of  this 
paper.  However,  the  rotational  curve  of  the  maser  disk  scales 
as  v  oc  r'031  (Greenhill  et  al.  1996),  so  even  a  factor  of  1.5 
decrease  in  the  circular  velocity  would  require  a  factor  of  ~4 
increase  in  the  radial  scale,  placing  the  CO-emitting  gas  well 
outside  of  the  maser  disk.  Alternatively,  it  would  be  possible  to 
place  the  HE  component  within  the  maser  disk  if  we  allow  that 
a  non-uniform  excitation  or  mass  distribution  enhances  the  CO 
emission  within  the  narrower  observed  velocity  range. 

The  strongest  argument  against  associating  the  HE  compo¬ 
nent  with  the  maser  disk  follows  from  a  comparison  of  the 
physical  parameters,  as  the  density  and/or  thermal  pressure  we 
estimate  for  the  HE  component  is  likely  to  be  smaller  than  that 
of  the  disk.  A  first  estimate  of  the  pressure  in  this  region  may 
be  obtained  by  considering  the  free-free-emitting  plasma  lying 
just  inside  the  maser  disk  ( r  ~  0.4  pc),  for  which  Gallimore 
et  al.  (2004)  estimate  neTe  ~  (6  x  105  cm  3)(6  x  106  K)  ~ 
1012  6  K  cm'3.  This  is  a  factor  of  ~103  larger  than  the  range 
P/kB  ~  10s  8— 109  8  K  cm'3  derived  from  our  LVG  modeling 
of  the  HE  component  (Table  3).  Separately,  Lodato  &  Bertin 
(2003)  have  derived  the  surface  density  profile  of  the  maser  disk 
required  to  reproduce  the  non-Keplerian  rotation  curve  and  find 
a  density  of  «h2  =  ( 1-5)  x  10s  cm'3  at  the  outer  edge  of  the  disk. 
Interferometric  mid-IR  observations  have  identified  a  structure 
of  hot  ( T  ~  800  K)  dust  with  a  size  of  ~0.45  x  1.35  pc,  which 
may  coincide  with  the  maser  disk  (Jaffe  et  al.  2004;  Raban  et  al. 
2009).  This  dust  temperature  sets  a  lower  limit  to  the  tempera¬ 
ture  of  the  concomitant  gas.  Densities  larger  than  108  cm'3  are 
in  excess  of  the  mh2  ~  105'9— 1071  cm'3  range  derived  from 
our  LVG  modeling  (Table  3),  and  models  simultaneously  re¬ 
quiring  77 H2  >  108  cm'3  and  7/ln  >  800  K  are  even  more  firmly 
excluded  (Figure  3). 

How  far  from  the  AGN  would  it  be  possible  to  find  molecular 
gas  at  the  relatively  low  pressure  we  attribute  to  the  HE 
component?  Neufeld  et  al.  (1994)  and  Neufeld  &  Maloney 
(1995)  have  investigated  the  properties  of  dense  gas  in  close 
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proximity  to  a  strong  X-ray  source  and  find  that  such  gas  is 
molecular  if  the  pressure  exceeds 

P/kB  >  I0n L43.5R~2N^-9  Kcm~3,  (3) 

where  1043  5L43.5  erg  s-1  is  the  2-10  keV  luminosity  and  the 
gas  is  Rpc  pc  from  the  source  and  is  shielded  by  a  column  of 
Nh  =  1024Af24  ctrT2.  For  NGC  1068,  we  estimate  L43.5  =  1 
(see  Section  7.2),  and  for  the  maser  disk  we  also  set  Rpc  =  1. 
Assuming  that  the  maser  disk  is  not  shielded  by  Compton-thick 
material  inward  of  1  pc  (i.e.,  setting  N24  ^  1),  this  expression 
supports  our  previous  conclusion  that  the  molecular  gas  in  the 
maser  disk  is  at  higher  pressure  than  our  HE  component.  If 
the  HE  emission  traces  the  Compton-thick  structure  blocking 
the  hard  X-rays,  then  by  construction  any  material  between 
this  gas  and  the  AGN  must  be  Compton-thin  (A24  <  1).  For 
P/kB  ~  10s  8— 109  8,  Equation  (3)  then  places  this  gas  at  a 
distance  of  Rpc  ~  4-13. 

We  conclude  that  the  HE  CO  emission  does  not  arise  from 
the  maser  disk,  primarily  due  to  the  large  difference  in  thermal 
pressures.  However,  the  observed  line  widths  and  the  results 
from  our  LVG  modeling  do  allow  us  to  construct  a  model  in 
which  the  HE  component  traces  the  nuclear  obscuring  material, 
provided  that  this  material  lies  at  least  a  few  parsecs  from 
the  AGN.  This  scale  may  be  broadly  consistent  with  clumpy 
torus  models,  in  which  the  obscuring  medium  is  composed  of 
clouds  distributed  from  the  sublimation  radius  (rsuh  <  1  pc)  to 
several  parsec  scales  (e.g.,  Honig  et  al.  2006;  Nenkova  et  al. 
2008).  These  models  are  currently  constrained  primarily  by 
IR  continuum  observations.  The  inclusion  of  a  gas  phase  in 
these  models  would  be  a  useful  next  step  to  further  evaluate  the 
prospects  for  attributing  the  FIR  CO  emission  to  gas  in  close 
proximity  of  the  AGN. 

7.2.  A  Comparison  with  the  Krolik  &  Lepp  (1989)  Torus  Model 

In  addition  to  the  detected  /upper  ±  30  transitions,  our  upper 
limits  to  the  /upper  ±  50  lines  provide  a  useful  constraint  on  any 
potential  high-excitation  nuclear  molecular  component.  Krolik 
&  Lepp  (1989)  modeled  the  molecular  emission  expected  from 
a  Compton-thick,  parsec-scale  torus.  They  calculated  the  CO 
emission  from  an  Ah  =  1024  cm-2  cloud  located  ~1  pc 
from  a  luminous  (Lx  ~  1044  erg  s_1)  hard  X-ray  source, 
in  which  the  heating  is  dominated  by  Compton  scattering  of 
10-100  keV  photons.  In  this  model  the  CO  SED  scales  as 
'/upper  UP  40  /upper  ~  58,  and  the  absolute  line  luminosities 
are  proportional  to  the  total  absorbed  10-100  keV  luminosity 
(/absL  v44;  in  units  of  1044  erg  s-1).  Several  of  our  upper  limits 
to  CO  transitions  with  /Upper  =  34-47  independently  place  an 
upper  limit  to  this  component  corresponding  to  /absC.v44  0.1, 
while  the  non-detection  of  the  CO(44-43)  line  achieves  the 
most  stringent  limit  of  /a bsCX44  <  0.09  (Figure  9(a)).  This  limit 
may  be  reduced  by  stacking  the  individual  non-detections.  We 
have  obtained  such  a  stack  by  first  scaling  the  spectrum  of  each 
undetected  line  by  (44//upper)2,  which  references  each  spectrum 
to  that  of  CO(44-43)  under  the  assumption  that  the  CO  fluxes 
scale  as  /3  .  We  then  calculated  a  weighted  average  of  eight 

lines  with  low  noise  that  fall  in  clean  spectral  regions  and  binned 
to  600  km  s“ 1 .  The  result  is  shown  in  Figure  1 1 .  The  upper  limit 
to  this  stacked  line  pushes  the  limit  of  the  Krolik  &  Lepp  (1989) 
model  to  /abs ^44  <  0-038  (Figure  9(a)). 

Iwasawa  et  al.  (1997)  and  Colbert  et  al.  (2002)  model  the 
reflected  X-ray  spectrum  of  NGC  1068  and  estimate  intrinsic 
luminosities  of  L2-iokeV  =  1043— 1043  7  erg  s  1  (corrected 


-4000  -2000  0  2000  4000 

Velocity  [km/s] 


Figure  11.  Baseline-subtracted  stack  of  eight  transitions  with  7Upper  =  34^47, 
for  comparison  with  the  Krolik  &  Lepp  (1989)  models.  Dashed  lines  mark  the 
drier  uncertainty. 

to  the  d  =  14.4  Mpc  used  here).  A  comparison  of  the 
measured  [O  in]  15007  line  luminosity  with  a  log([0  m \/Lx)  — 
—  1.89  ratio  established  for  type  1  Seyferts  yields  Z,2-iokeV  = 
1043'5  erg  s_1,  while  a  comparison  of  the  estimated  Lb0i  with 
Z-2-iokev/Cboi  ~  0.1  yields  L2-iokeV  ~  10436  erg  s'1  (Colbert 
et  al.  2002,  and  references  therein).  Melendez  et  al.  (2008) 
derive  a  relation  between  L[Oiv]26/im  and  Z-2-iokeV  f°r  a  sample 
of  hard  (14-195  keV)  X-ray-selected  type  1  Seyferts.  Using 
this  relation,  along  with  L[Oiv]26Mm  —  4.7  x  1041  erg  s'1 
(Sturm  et  al.  2002),  yields  L2-iokeV  ^  1044  erg  s-1.  Given 
this  range  of  estimates,  we  adopt  L2-tokeV  =  1043'5  erg  s-1  and 
increase  this  by  1.43  to  correct  to  the  10-100  keV  range  (for 
Lv  oc  v"1).  Assuming  that  25%  of  this  power  (for  a  covering 
fraction  /cov  =  0.25)  is  absorbed  in  the  phase  modeled  by  Krolik 
&  Lepp  (1989)  gives  an  expected  /a bsCx44  —  0.11,  a  modest 
factor  of  ~2.9  higher  than  our  measured  upper  limit.  Given  the 
uncertainties  in  the  AGN  luminosity  and  covering  factor  of  the 
Compton-thick  material  in  NGC  1068,  as  well  as  uncertainties 
in  the  model  itself,  this  factor  of  ~2.9  overprediction  is  too  small 
to  reject  the  model  of  a  parsec-scale  torus  envisioned  by  Krolik 
&  Lepp  (1989).  This  non-detection  may,  however,  provide  a 
useful  constraint  on  the  more  recent  set  of  detailed  clumpy  and 
dynamical  torus  models  (e.g.,  Honig  et  al.  2006;  Nenkova  et  al. 
2008;  Wada  et  al.  2009;  Schartmann  et  al.  2010). 

8.  SUMMARY  AND  FUTURE  WORK 

We  have  detected  11  CO  transitions  in  the  /upper  =  14-30 
range  from  the  central  10"  (700  pc)  of  NGC  1068  and  obtained 
sensitive  upper  limits  to  most  other  transitions  up  to  /upper  ±  50. 
These  are  the  first  extragalactic  detections  of  FIR  CO,  which 
represent  a  new  probe  of  excited  molecular  gas  in  Seyfert 
nuclei.  The  detected  transitions  are  modeled  as  arising  from 
two  different  components:  an  ME  component  close  to  the  galaxy 
systemic  velocity  and  an  HE  component  that  is  blueshifted  by 
~80  km  s-1.  Our  main  results  are  as  follows: 

1.  Using  a  two-component  LVG  model,  we  derive  nH2  ~ 
105'6  cm”3,  7kj„  ~  170  K,  and  MH2  ~  106,7  Me  for  the 
ME  component  and  «h2  ~  106'4  cm-3,  Tkl„  ~  570  K,  and 
/V/H 2  ~  10s  6  M0  for  the  HE  component.  The  la  uncer¬ 
tainties  on  the  derived  temperatures  are  ±(0.20-0.35)  dex, 
while  for  density  and  mass  this  is  ±(0.6-0. 9)  dex.  Ex¬ 
tending  the  measured  CO  line  SED  to  lower-  and  higher-/ 
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lines  would  likely  reduce  these  uncertainties,  as  would  a 
joint  analysis  with  H2,  OH,  H2O,  and  other  complemen¬ 
tary  molecular  tracers.  We  will  follow  both  approaches  in 
forthcoming  papers. 

2.  These  two  components  are  denser  than  the  gas  traced  with 
millimeter  CO  observations,  and  the  HE  (and  possibly 
the  ME)  component  is  also  warmer.  The  ME  component 
makes  a  non-negligible  contribution  to  the  nuclear  mass 
budget,  although  large  uncertainties  in  the  masses  estimated 
from  both  the  FIR  CO  lines  and  CO(  1-0)  prevent  a  more 
quantitative  statement. 

3.  Comparing  the  FIR  CO  line  profiles  with  those  of  high 
spatial  and  spectral  resolution  observations  of  CO(2-l)  and 
H2  1-0  S(  1 )  allows  a  first  estimate  of  the  origins  of  the  ME 
and  HE  components  within  the  central  10".  Good  matches 
are  found  with  H2  1-0  5(1),  which  for  the  ME  component 
suggests  an  origin  in  the  ~200  pc  diameter  H2  ring.  The 
blueshifted  HE  lines  may  also  be  consistent  with  an  origin 
in  the  bluest  regions  of  the  H2  ring  but  are  better  matched 
to  the  clump  of  infalling  molecular  gas  ~40  pc  north  of  the 
AGN. 

4.  The  ME  component  is  nicely  consistent  with  models  of 
X-ray  heating  of  gas  in  the  CND.  A  shock  model  is  also 
possible,  although  due  to  the  uncertainties  in  the  amount  of 
mechanical  power  available  for  dissipation  in  slow  shocks 
and  the  evidence  for  X-ray-driven  chemistry  in  the  CND, 
we  favor  the  X-ray-heating  scenario.  Far-UV  heating  is 
unlikely. 

5.  The  HE  component  is  also  consistent  with  either  X-ray-  or 
shock  heating.  X-ray  heating  would  best  fit  with  an  origin  in 
the  cloud  ~40  pc  north  of  the  AGN,  supporting  the  results 
of  the  line  profile  matching  (point  3).  Shocks  triggered  by 
the  interaction  of  the  radio  jet  with  this  clump,  or  arising 
from  the  H2  ring,  are  also  plausible.  Far-UV  heating  is 
unlikely. 

6.  The  thermal  pressure  of  our  HE  component  is  too  low  to 
be  attributed  to  gas  within  the  parsec-scale  H2O  maser 
disk  centered  on  the  AGN.  However,  the  pressure  may 
be  consistent  with  gas  located  ~4  pc  or  more  from  the 
AGN,  and  this  gas  could  potentially  provide  the  Ah  ~ 
1025  cm”2  column  obscuring  the  nuclear  hard  X-rays.  Our 
non-detections  of  /upper  =  34-47  lines  place  an  upper  limit 
to  the  Krolik  &  Lepp  (1989)  torus  model  that  is  a  factor 
of  ~2.9  lower  than  the  expected  signal,  although  due  to 
the  uncertainties  involved  in  applying  this  model,  this  non¬ 
detection  is  insufficient  to  rule  out  the  parsec-scale  torus 
paradigm.  The  inclusion  of  a  gas  phase  in  the  current  set 
of  clumpy  and  dynamic  torus  models  and  a  comparison  of 
the  predicted  CO  line  SED  with  our  detections  and  upper 
limits  would  be  a  useful  next  step. 
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